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Abstract 

R-trees can be used to store and query sets of point data in two or more dimensions. An 
easy way to construct and maintain R-trees for two-dimensional points, due to Kamel and 
Faloutsos, is to keep the points in the order in which they appear along the Hilbert curve. The 
R-tree will then store bounding boxes of points along contiguous sections of the curve, and the 
efficiency of the R-tree depends on the size of the bounding boxes—smaller is better. Since 
there are many different ways to generalize the Hilbert curve to higher dimensions, this raises 
the question which generalization results in the smallest bounding boxes. Familiar methods, 
such as the one by Butz, can result in curve sections whose bounding boxes are a factor fl(2 d ^ 2 ) 
larger than the volume traversed by that section of the curve. Most of the volume bounded by 
such bounding boxes would not contain any data points. In this paper we present a new way 
of generalizing Hilbert’s curve to higher dimensions, which results in much tighter bounding 
boxes: they have at most 4 times the volume of the part of the curve covered, independent of 
the number of dimensions. Moreover, we prove that a factor 4 is asymptotically optimal. 
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1 Introduction 


1.1 Space-filling curves and spatial index structures 

A d-dimensional space-filling curve is a continuous, surjective mapping from R. to R d . In the late 
19th century Peano m described such mappings for d = 2 and d = 3. Since then, various other 
space-filling curves have been found, and they have been applied in diverse areas such as spatial 
databases, load balancing in parallel computing, improving cache utilization in computations on 
large matrices, finite element methods, image compression, and combinatorial optimization mmm- 
In this paper we present new space-filling curves for d > 2 that have favourable properties for use 
in spatial data structures. 

In particular, we consider data structures for d-dimensional points such as R-trees [12) . In such 
data structures, data points are organised in blocks, often stored in external memory. Each block 
contains at most B points, for some parameter B , and each point is stored in exactly one block. 
For each block we maintain a bounding box, which is the smallest axis-aligned d-dimensional box 
that contains all points stored in the block. The bounding boxes of the blocks are stored in an 
index structure, which may often be kept in main memory. To find all points intersecting a given 
query window Q , we can now query the index structure for all bounding boxes that intersect Q; 
then we retrieve the corresponding blocks, and check the points in those blocks for answers to our 
query. We may also use the index structure to find the nearest neighbour to a query point q: if 
we search blocks in order of increasing distance from q , we will retrieve exactly the blocks whose 
bounding boxes intersect the largest empty sphere around q. The grouping of points into blocks 
determines what block bounding boxes are stored in the index structure, and in practice, retrieving 
these blocks is what determines the query response time jj]. 

If we store n points in d dimensions with B points in a block, 0((n/I3) 1-1 / d ) blocks may need to 
be visited in the worst case if the query window is a rectangular box with no points inside m, and 
Q(n/B) blocks may need to be visited if the query window is an empty sphere. The Priority-R-tree 
achieves these bounds )2j, whereas a heuristic solution by Kamel and Faloutsos da, which is 
explained below, may result in visiting Q(n/B) blocks even if the query window is a rectangular 
box with no points inside [2|. However, experimental results for (near-)point data and query ranges 
with few points inside [8] indicate that the approach by Kamel and Faloutsos seems to be more 
effective in practice for such settings. Moreover, regardless of the type of data and query ranges, a 
structure based on the ideas of Kamel and Faloutsos is much easier to build and maintain than a 
Priority-R-tree [2J. 

Kamel and Faloutsos proposed to determine the grouping of points into blocks as follows: we 
order the input points along a space-filling curve and then put each next group of B points together 
in a block (see Figure [ljb) ). Note that the number of blocks retrieved to answer a query is simply 
the number of bounding boxes intersected. Therefore it is important that the ordering induced by 
the space-filling curve makes us fill each block with points that lie close to each other and thus 
have a small bounding box. 

Kamel and Faloutsos proposed to use the Hilbert curve j9]j for this purpose. One way to describe 
the two-dimensional Hilbert curve is as a recursive construction that maps the unit interval [0,1] to 



Figure 1: (a) Sketch of Hilbert’s space-filling curve, (b) Blocks of an R-tree or similar data structure 
with B = 3. (c) Box-to-curve ratio of the section between p and q = area of the bounding box of 

the curve section S between p and g, divided by the area covered by S: 12 * 12/87 « 1.66. 
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the unit square [0, l] 2 . We subdivide the square into a grid of 2 x 2 square cells, and simultaneously 
subdivide the unit interval into four subintervals. Each subinterval is then matched to a cell; thus 
Hilbert’s curve traverses the cells one by one in a particular order. The mapping from unit interval 
to unit square is refined by applying the procedure recursively to each subinterval-cell pair, so 
that within each cell, the curve makes a similar traversal. The traversals within these cells are 
rotated and/or reflected so that the traversal remains continuous from one cell to another (see 
Figure Qa)). The result is a fully-specified mapping / : [0,1] —> [0, l] 2 from the unit interval to 
the unit square. The mapping is easily reversed, and thanks to the fact that the curve is based 
on recursive subdivision of a square in quadrants, the reversed mapping can be implemented very 
efficiently with coordinates represented as binary numbers. This gives us a way to decide which of 
any two points in the unit square is the first along the curve. 

We can sketch the shape of the curve by drawing, for the fc-th level of recursion, a polygonal 
curve, an approximating curve Ak, that connects the centres of the 4 fc squares in the order in which 
they are visited. In fact, the mapping / can also be described as the limit of the approximating 
curves Ak as k goes to infinity. Explicit descriptions of the approximating curves help us to reason 
about the shapes of curve sections, and thus, about the extents of their bounding boxes. For ease 
of notation, in this paper we scale the approximating curve for any level A; by a factor 2 k and 
translate it so that its vertices are exactly the points {0, - -., 2 fc — l} 2 . 

A d-dimensional version of Hilbert’s curve could now be described by a series of curves Af. 
for increasing k, each visiting the points {0,..., 2 k — l} d , where each point corresponds to a 
d-dimensional cube of width l/2 k in the unit hypercube. For d > 3, there are many ways to define 
such a series of curves Huai, but their distinctive properties and their differences in suitability 
for our purposes are largely unexplored. 

1.2 Our results 

In this paper we present a family of space-filling curves, for any number of dimensions d > 3, with 
two properties which we call well-foldedness and hyperorthogonality —Hilbert’s two-dimensional 
curve also has these properties. We show that these properties imply that the curves have good 
bounding-box quality as defined by Haverkort and Van Walderveen |7j. 

More precisely, for any 0 < a < b < 1, let /([a, 6]) denote the section of the space-filling curve 
/ from /(a) to /(d), that is, /([a, b}) = U a <t<&{/(*)}- 

The box-to-curve ratio (BCR) of a section f([a,b]) (denoted BCR(/([a, &]))) is the volume of the 
minimum axis-aligned bounding box of /([a, b\) divided by the volume (d-dimensional Lebesgue 
measure) of /([a, b]), see Figure[ljc). 

The worst-case BCR of a space-filling curve / is the maximum BCR over all sections of /. We 
show that the worst-case BCR of a well-folded, liyperorthogonal space-filling curve is at most 4, 
independent of the number of dimensions. Moreover, we show that this is asymptotically optimal: 
we prove that any d-dimensional space-filling curve that is described by a series of curves Ak as 
defined above, has a section with BCR at least 4 — 0(l/2 d ). In contrast, the d-dimensional “Hilbert” 
curves of Butz [4], as implemented by Moore m , have sections with BCR in fl(2 d / 2 ). 

In Section [l.3| we introduce basic nomenclature and notation. Section [2] defines the concept 
of well-foldedness, and presents sufficient and necessary conditions for approximating curves of 
well-folded space-filling curves. Section [3] introduces the concept of hyperorthogonality. We present 
sufficient and necessary conditions for approximating curves of well-folded space-filling curves to 
be hyperorthogonal. The necessity of these conditions is then used to prove that any section of a 
hyperorthogonal well-folded space-filling curve has good box-to-curve ratio. Our next task is to 
show that hyperorthogonal well-folded curves actually exist, and this is the topic of Section [4] 

We combine the conditions from the previous sections to learn more about the shape of 
hyperorthogonal well-folded curves, and in particular about self-similar curves (Section [5j. It turns 
out that in two, three, and four dimensions, there are actually very few self-similar, well-folded, 
hyperorthogonal curves (Corollary |48| ); in five and more dimensions, more such curves exist. In 
Section [bj we make a few remarks about how to implement a comparison operator based on 
self-similar, well-folded, hyperorthogonal curves in any number of dimensions greater than two. 
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Finally, in Section [7] we compare the bounding box quality of hyperorthogonal well-folded curves 
to lower bounds and to the bounding box quality of Butz’s generalization of Hilbert curves, and 
we discuss directions for further research. Pseudocode for the comparison operator discussed in 
Section [6] is given and explained in Appendix [A} 

1.3 Nomenclature and notation 

General notation 

— By D we denote 2 d . 

— By sign(i) we denote the sign of i, that is, sign(z) = — 1 if i < 0; sign(i) = 0 if i = 0, and 
sign(z) = 1 if i > 0. 

— By isneg(z) we denote the function defined by isneg(z) = 1 if i < 0, and isneg(z) = 0 if i > 0. 
Notice sign(z) = 1 — 2 * isneg(i). 

Vertices, edges, directions and axes 

— The universe in this article is the integer grid in d dimensions Z d . 

— A vertex is a point v = (u[l], v[2 ],... ,u[d]) £ Z d . 

— An edge e is an ordered pair of vertices (v,w) with distance |\w — u| | = 1. 

— The direction of an edge e = ( v , w ) is the number i £ {— d, —d + 1,..., d — 1, d} \ {0} such that 
w;[|z|] - n[|z|] = sign(z) and w[j] = v\j] if j ^ |z|. 

— The axis of an edge is the absolute value of its direction. Note that the edges (v, w) and (w, v) 
have opposite directions, but the same axis. 

— In our figures we will use horizontal lines for axis 1, vertical lines for axis 2 and lines with another 
orientation for axis 3. 

— By (ei, 62 , ■ • ■) we denote a path of edges with directions ei, e 2 ,.... 

Curves, length, volume, entry and exit 

— For the purposes of this paper, a curve is a curve on the grid , which is an ordered set of unique 
vertices where each subsequent pair of vertices forms an edge as defined above. Note that a curve 
never visits the same vertex more than once. Since a vertex and a direction determine an edge, a 
curve can alternatively be specified by the starting point and the listing of the directions of its 
edges in order. Note that curves are directed. 

— A space-filling curve is always a mapping / : [0,1] —>• [0, l] d , while any other curve discussed in 
this paper will be assumed to be a curve on the grid. 

— A free curve is a curve without a starting point, so with unspecified location: it is described by 
the directions of its edges only. 

— The reverse C of a free curve C is obtained by reversing the order of the edge directions and 
reversing the directions themselves, which means negating them. 

— The length of a curve is the number of edges, the volume of a subset of the grid is the number of 
vertices it contains. So the volume vol(C) of a curve C is its length + 1. 

— The first vertex of a curve is called the entry, the last vertex is called the exit. 

/.-Curves and k- cubes 

— A k-cube is a d-dimensional cube with 2 d * k points, so with a side of length 2 k — 1. 

— A k-curve is a Hamiltonian path on the integer grid in a k-cube. 

— Since each of the (integer) points of the cube is visited by the curve exactly once, its volume is 
2 d * k and the length of a /c-curve is 2 d * k — 1. 

Approximating curves 

— The space-filling curves under study in this paper will be approximated by curves on the grid as 
just defined. By Aq, A±,... we will denote a sequence of curves that approximates a d-dimensional 
space-filling curve, where Aq is a single vertex and Ak is a fc-curve. 

— By t’fcq, Vk. 2 , ■ ■ ■ ,Vk,K, where K = 2 d * k , we denote the vertices of Ak in order, and by ek,i we 
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Figure 2: Left: A parent curve A\ is inflated to create A 2 , which is composed of the child curves 
Cl,!, (71,2, Ci, 3 and (7i,4, and edges of A\ which are translated such that they connect the child 
curves to each other at their end points. Right: G( 3) with the directions of its edges. 


denote the direction of the edge (vk,i,Vk,i+ 1 ). 

— Each vertex Vk,i of Ak represents a d-dimensional hypercube Hk,i of width l/2 fc that is visited 
by the space-filling curve approximated by Ak- The vertices Vk+i,D*i-D+i, ■ ■ ■ ,14+1 ,D*i of Ak+ 1 
model the order in which the space-filling curve traverses the d-dimensional hypercubes of width 
l/ 2 fc+1 whose union is Hki- 

— Therefore it must be possible to construct Ak +1 from Ak, which we call the parent curve, by 
inflation: we replace each vertex 14 ,, of the parent curve with a 1 -curve Ck,i (a child curve), whose 
vertices are those of the unit cube, translated by 2 * 14 , 3 . Each edge (vk,i, 14 , 3 + 1 ) of the parent 
curve is replaced by an edge (14+1,D*i, 14+1,o*i+i) in Ak +1 of the same direction, connecting the 
exit of Cfc,j to the entry of Ck,i+ 1 , see Figure [2j left. Note that not just any choice of child curves 
results in a valid (k + l)-curve. The 1-curves that replace the vertices have to be chosen carefully 
such that for each edge ( 14 ,*, 14 , 1 + 1 ) of the parent curve, there is indeed an edge in the grid from 
the exit of Ck,i to the entry of Ck,i+ 1 - In Section 2.3 we will discuss how the 1-curves should be 
constructed so that they match up. 

In what follows, the first subscripts to v, e, H and C will usually be omitted if they are clear 
from the context, for example, if k is fixed, or if the approximating curve in question is otherwise 
specified. For example, “a child curve C, of Ak- 1 ” should be read as: “a child curve Ck- i,»” (and 
it would be a subcurve of Ak). 

Observe that our definition of curves on the grid restricts the generalizations of Hilbert curves 
under study to face-continuous curves, that is, each pair of consecutive d-dimensional hypercubes 
along the space-filling curve must share a (d — 1)-dimensional face. In Section 7.2 we will discuss 
why, in the context of this paper, this restriction is justified. 


2 Well-folded curves 

2.1 Gray codes and definition of well-folded curves 

In the process of inflating, we will restrict ourselves in this paper to replacing vertices with isometric 
images (like translations, rotations, reflections or taking the reverse, shortly: all distance-preserving 
mappings) of one particular 1-curve, namely the free curve G(d) that follows the so-called binary 
reflected Gray code. 

Definition 1. The free curve G(d) is defined recursively as follows: G{ 0) is empty; G(d) is the 
concatenation of G[d — 1), (d), and G(d — 1). 

For example, G(2) is the free curve (1, 2, —1) (Figure [2j left), G( 3) is shown in Figure [2] right, 
and G( 4) is the free curve (1, 2, —1, 3,1, —2, —1,4,1,2, —1, —3,1, —2, —1). 
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The length of G(d) is, by induction, 2 d — 1, which is the maximum length of a Hamiltonian 
path on the unit cube in r L d . 

The following properties of G(d) are well-known: 

Lemma 2. If d> 2, then, in G(d) as well as G(d), edges with axis 1 and edges with other axes 
alternate, starting with (1) and ending with (—1). 

Proof. Straightforward by induction on increasing d , with base case d = 2. □ 

Lemma 3. Let G°(d) be G(d) with entry point (0,..., 0), where d > 1. Then the vertices of G°(d) 
are those of the 1-cube {0, l} d and the exit point is the point v = (0,..., 0,1) with v[d] = 1 and 
v\j] = 0 for j < d. 

Proof. We can prove this by induction on increasing d. The base case d = 1 is easy to verify. 
Now suppose the lemma holds for d — 1, that is, the vertices of G°(d — 1) are those of the 1-cube 
{0, l}* 1 ” 1 and the exit point is (0,..., 0,1). Or to put it differently, the exit point of G°(d — 1) is 
the endpoint of an edge (d — 1) that starts in the origin. 

Now recall that G(d) is the concatenation of G(d— 1), (d) and G{d — 1). Therefore the exit 
point of G°(d) can be found as the endpoint of a curve ((d — 1), d, — (d — 1)), starting in the origin. 
Clearly, this endpoint is (0,..., 0,1). The vertices of G°(d) are those of G(d — 1) starting in the 
origin plus those of G(d — 1) ending at (0,..., 0,1), or equivalently, those of G(d— 1) starting in the 
origin plus those of G(d — 1) starting at (0,..., 0,1). Together these constitute the set {0, T} d . □ 

The following lemma will prove useful in Sections [3] and |7.1| but can be skipped on first reading: 

Lemma 4. The axes of the first (and last) n edges of G{d) constitute the set {1,... , to}, where 
m = 1 + |k>g 2 (n)J = flog 2 (n+ l)]. 

Proof. For 1 + [log 2 (n)J = to = [log 2 (n + 1)] we have 2 n > 2 m > n + 1 and thus, n > 2 m ~ 1 and 
n < 2 m — 1. Therefore the first n edges of G(d) include at least a full G(m — 1) and an edge (to), 
and not more than a full G(m). In a symmetric way, the last n edges of G(d) include at least an 
edge (—to) and a full G(?n — 1), and not more than a full G(to). It follows that the first or last n 
edges cover to different axes. □ 

Definition 5. A curve is well-folded if it is a single vertex, or if it is obtained by inflating a 
well-folded curve by replacing its vertices by isometric images ofG(d). A space-filling curve is 
well-folded if its approximating curves are well-folded. 

Note that in two dimensions, all possible 1-curves are in fact isometric images of G(2), so any 
face-continuous space-filling curve based on recursive subdivision of a square into four squares must 
be well-folded (for example, Hilbert’s curve or the /3H-curve [17]). 

In higher dimensions, the most common generalizations of the Hilbert curve are well-folded as 
well, but there are also face-continuous curves based on recursive subdivision of a cube into eight 
cubes that are not well-folded (using generators of types B and C from Alber and Niedermeier [1:5]). 
In Section]?] we will briefly get back to non-well-folded curves; until then, we will focus on well-folded 
curves. 

2.2 Notation for isometries of Gray codes in well-folded curves 

The isometric transformations of 1-curves which we need in this paper are those of the hyperocta- 
hedral group of symmetries of the hypercube. This group is the product of the symmetric group 
Sd (the group of all permutations of the d coordinate axes) and the group of 2 d reflections formed 
by all combinations of reflections in hyperspaces orthogonal to the coordinate axes. Thus there are 
d\ * 2 d such transformations. 

To distinguish these transformations, we will use signed permutations. A signed permutation n 
is a bijection from {— d, —d + 1,..., d — 1, d} \ {0} to itself with the property that n(—k) = —i r(fc) 
for Jc e {1,..., d}. It is denoted by [7r(l), 7t(2), ..., 7r(d)]. 
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Given a fc-cube H , a signed permutation 7 r specifies the isometry that maps H onto itself and 
maps the direction k to the direction n(k). If tt is a signed permutation, then 7 r(A’) denotes the 
application of 7 r to all elements of the vector, set, or sequence A; | 7 r| denotes the permutation 
[|7r(l)I, |7 t( 2)|, ..., |7r(d)|]; and 7r _1 denotes the inverse of 7r, that is, 7r _ 1 (a;) = y if and only if 
n(y) = x. 

Note that signed permutations do not allow us to express the isometric transformation that 
consists of reversing a curve. For now, this is not a problem, because the reversal of G(d) is 
identical to its reflection in coordinate d. 

We define the orientation of an isometry of G(d) as the direction of the vector from entry to 
exit. A direct corollary of Lemma [3] is the following: 

Corollary 6 . Let G°(d) be G(d) with entry point (0,... ,0), and let n be a signed permutation. 
Then the coordinates of the entry point a of n(G(d)) are given by a[j] = isneg (tt - 1 (j)) for 
j £ {1,..., d}; the orientation of Tr(G(d)) is n(d); and the coordinates of the exit point b of i:(G{d)) 
satisfy b[j] = 1 - a[j] for j = \n(d)\ and b[j\ = a[j\ for j 7 ^ |tt(o?)|. 

Isometries in approximating curves Consider a sequence of well-folded approximating curves 
Aq,Ai,.... By ak,i we denote the transformation (modulo translation) that is applied to G(d) to 
obtain the 1-curve 6 \-p that replaces vertex Vi of Aj- in the inflation of A/ c to Ak+ For example, 
for the curves in Figure[2j left, we have (Tip = [—1, 2]; cri i2 = [—2,1]; cry 3 = crip = [2, —1]. As with 
v, e, H and C, the first subscript will usually be omitted if it is clear from the context. 

2.3 Conditions on edges and isometries in well-folded curves 

As observed before, when inflating a curve A*,, the 1-curves that replace the vertices of Ak have to 
be chosen carefully such that the exit of C, and the entry of Ci +1 constitute an edge with direction 
ei. For this we need the conditions as stated in Theorem [7] below. 

Theorem 7. Given a well-folded approximating curve A & for a particular, fixed level k > 0. 
Inflating Ak to Ak +1 results in a well-folded approximating curve Ak+i if and only if, for each 
1 <i< 2 d * k : 

• for j £ { 1 ,..., d} we have sign (crF^j)) = sign if and only if j equals neither or 

both of \cri(d)\ and |e;|; otherwise sign (cr^j)) = - sign (af 1 ^)); 

• sign (a-^id)) = 1 . 

Proof. By construction, Ak+\ is obtained by replacing the vertices of Ak by isometric images of 
G(d). The challenge is to prove that the above conditions are necessary and sufficient to guarantee 
that Ak +1 is indeed a curve, and hence, well-folded. 

Recall that C z = a-i(G(d)) is the 1-curve (modulo translation) that replaces vertex Vi of Ak in 
the process of inflation and = (vi,Vi + 1 ). For a given i, let a , b and c be, respectively, the entry 
of Ci, the exit of Cj, and the entry of Ca+ 1 , all relative to the point 2 * Vi- By Corollary [ 6 ] we 
have a[j] = isneg (cr ~ (j)) and c[j] = isneg (<JT_i(j)) (mod 2) for j £ {1,... ,d}. Our task is to 
establish the conditions under which (b,c) is indeed an edge with direction e,, that is: for j = |e;| 
we should have c[j] = b[j] + sign(ei), and for j 7 ^ |e,;| we should have c[j ] = b[j]. 

Note that if ( 6 ,c) is an edge with direction ei, then the path ( <Ji{d),ei) brings us from a 
via b to c. Since each edge increments or decrements one coordinate by 1, it follows that if 
and only if j equals neither of both of \ai(d)\ and \ef\, we have c[j] = a[j] (mod 2 ) and thus, 
isneg (crf 1 ^)) = isneg (crf+iij)). This proves that the first condition of the theorem is necessary. 

Conversely, the first condition of the theorem, together with Corollary[ 6 j gives us that b[j] = c[j ] 
(mod 2) if and only if j 7 ^ |e» |, as witnessed by the following table. The third column is equivalent 
with the first condition of the theorem, as derived in the previous paragraph. The fourth column 
follows from Corollary [ 6 ] The fifth column follows from the third and the fourth. 
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j = \ a i(d)\ 

j = M 

a[j] = c \j\ ( mod 2 ) 

a[j] = b[j] (mod 2) 

b\j] = c[j ] (mod 2) 

no 

no 

yes 

yes 

yes 

no 

yes 

no 

yes 

no 

yes 

no 

no 

no 

yes 

yes 

yes 

yes 

no 

no 


In fact, since i>i[j] and n,;+i[j] are equal if j ^ \ei\, we get b\j] = c[j] (without (mod 2)) if j ^ \a\. 
On the other hand, if j = |e» |, we have, so far, only established b[j] ^ c[j] (mod 2). 

To complete the proof of Theorem [?J we will now show that, given j = |e*| and c[j\ ^ b[j] 
(mod 2), we actually have c[j] = b[j] + sign(ei) if and only if the second condition of the theorem 
is satisfied. In fact, the second condition, sign (<rf\ (ej) = 1, expresses that, within the 1-cube 
filled by Cj+i, the entry is on the side that is adjacent in direction e, to the 1-cube filled by Ci. To 
analyse this in more detail, we distinguish two cases: first, sign(e,) = 1, and second, sign(ei) = — 1. 

If sign(ei) = 1, then Vi + \[j] = Vi[j] + 1. Hence, by the fact that Ak+\ is obtained by inflation 
from Ak , we have 0 < b[j] < 1 and 2 < c[j] < 3. Moreover, given c[j] b[j] (mod 2), we have 
that (6[j],c[j]) is either (0,3) or (1,2). So we have c[j] = b[j] + sign(ei) = b[j] + 1 if and only 
if c[j] = 2 = 0 (mod 2). By Corollary this is the case if and only if isneg (cr,^ (j)) = 0; with 
j = \et\ = et we can rewrite this as sign = 1. 

Similarly, if sign(ej) = —1, then u*+ib1 = Vi[j] — 1. Hence, by the fact that Ak+i is obtained by 
inflation from Ak , we have 0 < b[j] < 1 and —2 < c[j] < — 1. Moreover, given c[j] ^ b[j] (mod 2), 
we have that (b[j\, c[j]) is either (0, —1) or (1, —2). So we have c[j] = b[j] + sign(ej) = b[j] — 1 if and 
only if c[j] = —1 = 1 (mod 2). By Corollary 6 this is the case if and only if isneg (crC 1 -, (j)) = 1; with 
j = |ei| = -e, we can rewrite this as sign = sign (er”^—j)) = - sign (eq+^O')) =1. □ 

Given the edges of Ak and the signs of the inverse permutations, Theorem [T] allows us to 
determine the last elements (Ji{d) of each permutation. Conversely, given the edges and the last 
elements of each permutation, Theorem [7] allows us to determine the signs of each permutation. Note 
that this leaves d — 1 elements of each |cr.j | unspecified and without consequence: any permutation 
of those elements will do. 

Observation 8. Let f be a well-folded space-filling curve approximated by Aq, A±,..and let 
x = /(0) be the starting point of f. Then x[j] = i sne g( CT fc i(j))/2 fe+1 - 

In other words, the digits of the binary representation of x[j] behind the fractional point are 
isneg (<7(7 1 (j )), isneg (af J (j)), isneg {af {(j)),.... 

3 Hyperorthogonal well-folded curves 

So far, we have been defining and discussing properties of curves that are in fact common to 
the previously best-known generalizations of Hilbert’s curve to higher dimensions. We will now 
introduce a new property that is not satisfied by any of the previously known generalizations that 
we are aware of, and which will prove useful in designing novel curves with good box-to-curve 
ratios. 

3.1 Definition and characterization 

Definition 9. We call a curve hyperorthogonal if and only if, for any n £ {0, ...,d — 2}, 
each sequence of 2 n consecutive edges have exactly n + 1 different axes. A space-filling curve is 
hyperorthogonal if its approximating curves are hyperorthogonal. 

Notice that an n-dimensional 1-cube (in Z d ) can hold at most 2™ — 1 consecutive edges of a 
curve, so any curve constructed by inflation contains sets of 2” edges that have at least n + 1 
different axes, for each n < d — 1. Hyperorthogonality requires that this holds for every set of 
2" edges, provided n < d — 2. The definition leaves little room for being made more strict (see 
Inset [I]). 
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Inset 1 No room for a stricter definition of hyperorthogonality 

In Definition [9j the upper bound on n cannot be raised to d — 1, as this would require that, in two 
dimensions, no pair of consecutive edges would have the same direction. It is easy to see that a 
2-curve with this property cannot be constructed. 

Consider a square of four by four vertices, that is, a 2-cube for d = 2. Within this square, 
four vertices lie in a corner. Of these four corner vertices, let s be the second one visited by the 
curve. The curve must visit at least two vertices before s aud at least two vertices after s in order 
to reach the other corners of the square. Let S be the unit square (quadrant of four vertices) that 
contains s. 

Now consider the sequence that consists of the two edges that precede s and the two edges 
that follow s. Since this sequence visits five vertices, it clearly does not fit in S, and therefore the 
2 d ~ 1 = 2 edges preceding s or the 2 d ~ 1 = 2 edges following s do not fit in a two-dimensional unit 
cube. Hence, either two edges preceding s or the two edges following s must be collinear. 

Hyperorthogonality still allows that less than 2 n , but more than 2 n ~ 1 , consecutive edges also 
span a (n + l)-dimensional space—this is also necessary, since otherwise, even if d > 3, any three 
consecutive edges would be restricted to alternating between two dimensions, which, by induction, 
would restrict the whole curve to edges alternating between two dimensions. 


For d = 2, hyperorthogonality requires only that each single edge spans a one-dimensional 
space, which is obvious. So all two-dimensional curves are hyperorthogonal. 

For d = 3 each two consecutive edges must span a two-dimensional space, so each pair of 
consecutive edges must be orthogonal. (For that reason the property is called ‘hyperorthogonal’ 
for higher dimensions as well.) 

Note that G(d) is hyperorthogonal for all d. 

As can be seen by inspecting familiar generalizations of Hilbert curves to three dimensions, if we 
construct a sequence of curves Aq, ..., A k in three or more dimensions by inflation, using isometric 
images of G(d) to inflate vertices, then Ak is not necessarily hyperorthogonal, even though G(d) is 
(see, for example, the Butz-Moore curve in Figure [5j right, where there are two collinear edges 
along the top back edge of the cube). The next theorem states what conditions the isometries 
should fulfill in order to obtain hyperorthogonal curves. 

Definition 10. The depth of a direction a in a signed permutation 7 r, denoted depth(7r, a), is 
defined as follows: if |a| £ {|7r(d)|, |7r(d — 1)|}, then depth(7r, a) = 0, otherwise depth(7r,a) is the 
number j such that \tt( d— 1 — j)\ = |a|. 

So the depth of n(d — 2) is 1, the depth of 7r(l) is d — 2, and since each axis occurs in 7r, each 
direction has a depth in 7 r. 

Theorem 11. For fixed k, let K = 2 d * k , and let Aq, ..., A k +i be a sequence of well-folded curves 
constructed by inflation (with all the associated notation introduced in the previous sections). 
Suppose Ak is hyperorthogonal , then A k +\ is hyperorthogonal as well if and only if the following 
conditions are satisfied: 

1. for each i £ {1,..., K - 1}: depth(cr fcji , e k ,i) = 0 = depth((T fe4+ i, e k ,i); 

2. for each i £ {1,..., K — 1} and each direction a: 

| depth(cr fcj i, a) - depth(a- fcji+ i, a)| < 1. 

Proof. As usual, we will omit the subscripts k in this proof. 

Necessity: Suppose condition [l] is violated, that is, |e^| £ {|cq(l)|,..., |afid — 2)|} or 
N e {|«T i+1 (i)|,...,| <Ji+i(d — 2)}. We analyse the first case \ef £ {|cq(l)|,..., \afid — 2)|}, the 
second case is symmetric. 

Consider the last 2 d_2 —1 edges of afiG^d)): by Lemma[4| their axes form the set {|cr.j(l)|,..., |cq(d— 
2)|}. In Afc + i, these edges will be followed by an edge with axis |e,;| £ {|cq(l)|,..., |afid — 2)|}. 
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Thus we get a sequence of 2 d ~ 2 edges with only d — 2 different axes: too few for hyperorthogonality 
as defined by Definition [9] 

Now suppose condition |T] is satisfied, but condition [2] is violated, that is, there are h and j such 
that |crj(h)| = \(Ji + \{j)\ or |<r i (j)| = \a i+ i{h)\ and h+l<j<d — l or h + 2 < j = d. We analyse 
the first case |er,;(/i)| = |eq + i(j)|, the second case is symmetric. 

Observe that |<jj(l)|,..., \<Ji(h)\ all differ from \a\, since h < d — 3. Also, in <7j+i(G(d)), the axes 
|(Tj + i(l)|,..., \cr i+ i(h + 1)| all differ from |e»| as well as from |cr i+ i(j)| = |<7j(/i)|. Now consider the 
sequence of 2 h+1 edges that consists of the last 2 h — 1 edges of eq(G(d)), followed by the edge with 
direction e,;, and the first 2 h edges of cr i+ i(G(d)). By Lemma |ij the last 2 h edges of this sequence 
have h + 1 different axes |o"i+i(l)|,..., \a i+ i(h, + 1)|, while the first 2 h edges contribute two more 
different axes, namely \<Ji(h)\ and |e»|. Thus there are h + 3 different axes in this sequence of 2 h+1 
edges: too many for hyperorthogonality as defined by Definition [9] 

Sufficiency: We distinguish two cases: a sequence of 2 n edges E, with n < d — 2, either lies 
within a single 1-curve <7j(G(d)), or not. 

In the first case, let j be the highest index such that E includes an edge with axis <Ji(j). From 
Definition [l] we get that the axes of the edges of E are, in order, for some m < 2 n — 1, those of 
the last m edges of cr,;(G(j — 1)), followed by cr,;(j) and the axes of the first 2 n — m — 1 edges of 

tJi(G(j — 1)). Since either m or 2" — m — 1 must be at least 2 n ~ 1 , it follows from Lemma k] that 
the axes of these edges are exactly {|cx^(1)|,..., |cr,;(n)|, |cr i (j)|}, since j > n. Hence, the edges of E 
have exactly n + 1 different axes and satisfy the conditions for hyperorthogonality. 

In the second case, E consists of the last m edges in a 1-curve Uj(G(d)), followed by (e*), and 
2 n — ni — 1 edges in cq + i(G(d)). Assume m > 2 n — m — 1 (the opposite case is symmetric), and 
hence, 2 n ~ 1 < m < 2 n — 1 and 2™ — m — 1 < 2 n_1 — 1. By Lemma[4j the first in edges have axes 
{|er,(l)|,..., |<7j(n)|}. Since \ei\ £ {|c i(d — 1)|, |cq(d)|} and d — 1 > n, the edge (e,) contributes one 
more axis. The remaining edges have axes from {|cr 2 : + i(l)|,..., (n — 1)|}, which, because of 
the second condition of the theorem, is a subset of {|cr*(1 )|,..., |cq(ro)|}; hence these edges do not 
contribute any more axes. In total, the edges of E have exactly n + 1 different axes and satisfy the 
conditions for hyperorthogonality. □ 

3.2 Box-to-curve ratio < 4 

To bound the box-to-curve ratio (bcr) of sections of hyperorthogonal well-folded space-filling 
curves, we will make use of the following lemma: 

Lemma 12. For any n £ {0,1,..., d — 2}, each sequence of 2 n consecutive edges of a well-folded, 
hyperorthogonal curve lies inside an {n + 1)-dimensional unit cube. 

Proof. Definition [9] states that each sequence E of 2” consecutive edges of a hyperorthogonal curve 
lies inside an axis-aligned box that has non-zero width in exactly n + 1 dimensions. Therefore, to 
prove the lemma, we only have to show that in none of these n + 1 dimensions, the width is more 
than 1. 

We can prove this by inspection of the sufficiency proof of Theorem El In the first case, the 
width is not more than 1 in any dimension, since all of E lies inside a single unit cube. In the 
second case, E lies in the union of two unit cubes, which is a box with width 3 in dimension |e* |, 
and width 1 in the remaining dimensions. However, as the proof argues, E contains only one edge 
with axis |e* |; hence the width in this dimension is only one. □ 

Theorem 13. The box-to-curve ratio of any section of a hyperorthogonal well-folded space-filling 
curve is at most f. 

Proof. Consider a section s of a hyperorthogonal well-folded space-filling curve /, approximated by 
a series of curves A 0 , Ai,.... Let E^ be the subcurve of A *. that contains all vertices trepresenting 
hypercubes Hi of width l/2 fc , whose interiors are intersected by s. For {vh ,..., Vj} = Ek, the 
bounding box of s is contained in the smallest axis-aligned box that fully contains all hypercubes 
H h ,... ,Hj. 
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Now let k be the smallest index such that Ek +1 contains at least one vertex that represents a 
hypercube of width l/2 fe+1 that is fully contained in s. By this choice of k, the subcurve Ef. of A 
does not contain any vertex that represents a hypercube of width l/2 k that is fully contained in s. 
Thus, Ek contains only a single vertex x, or two vertices x and y. and E^+i consists of vertices 
from the respective child curves C x and C y that replace x and y in the inflation from A & to A^+i- 
Note that this implies that the bounding box of E^+i has at most the volume of two 1-cubes, 
that is 2 d+1 . Define E = E^+ 1 , let A be the maximum common subcurve of C x and E, and, 
if y exists, let Y be the maximum common subcurve of C y and E, otherwise Y = 0. Thus, 
vol(y) = vol(-E) — vol(A); without loss of generality, assume vol(Y) < vol(A'). Furthermore, let 
c = \e m i n ( Xt y) | be the axis of the connecting edge of A' and Y. 

A number of cases with smartly chosen boundaries for vol (E), vol(X) and vol(Y) can now be 
distinguished, as shown in the table below. In each case, we derive an upper bound MaxBoxVol 
on the bounding box volume, and a lower bound MinCrvVol on the number of vertices of E that 
represent hypercubes completely covered by s (this is usually all of E except for the first and last 
vertex). From this we can derive that the box-to-curve ratio is less than MaxBoxVol/MinCrvVol < 
4. 


Case 

MaxBoxVol 

MinCrvVol 

A 

2 d ~ 1 + 2 < vol {E) < 2 d+1 

2 d +i 

2 d ~i 

B 

2^-2 + 2 < vol(£) < 2 d ~ x + 1 and... 



Bl 

...and vol(Y) < vol(A) < 2 d ~ 2 

2 d 

2 d ~2 

B2 

...and 2 d ~ 3 < vol(Y) < 2 d ~ 2 < vol(A) 

3 9 d 

2 * Z 

3 r\d—2 

2 ' Z 

B3 

...and 1 < vol(Y) < 2 d ~ 3 

2 d 

2 d ~2 

B4 

...and vol(Y) = 0 

2 d 

to 

0- 

1 

to 

C 

3 < vol (E) < 2 d ~ 2 + 1 

4(vol(£0 - 2) 

vo1(f;) - 2 

D 

vol(-E) < 2 

2 

1 


Note that Bl, B2, B3, and B4 are subcases for the same bounds on vol(Fl), where B1 is the case of 
having small X, and B2, B3, and B4 are the cases of large X with various bounds on the size of Y. 
For cases A, B4, and D the bounds on the bounding box volume are trivial; cases Bl, B2, B3, 
and C require a more careful analysis. 

Bl: By Theorem[lI] for the axis c of the connecting edge between A' and Y we have depth(cr x , c) = 
0. Since vol(A) < 2 d ~ 2 , and thus, the length of X is at most 2 d ~ 2 — 1, Lemma[4]now tells us 
that the edges of A' have axes from |<t x (1)|, ..., \a x (d — 2)|, hence not including c. Therefore 
X is included in the half-cube ((d — l)-dimensional 1-cube) that consists of the vertices of 
C x that are adjacent to vertices of C y . Likewise, Y is included in the half-cube that consists 
of the vertices of C y that are adjacent to C x . These two half-cubes together constitute a 
d-dimensional unit cube of volume 2 d . 

B2: As in case Bl, Y is included in the half-cube that consists of the vertices of C y that are adjacent 
to C x . This half-cube, together with C Xl has a bounding box of volume | ■ 2 d . The minimum 
curve volume MinCrvVol is at least vol (E) — 2 = vol(A) +vol(F) — 2 > 2 d ~ 2 + 2 d ~ 3 = | • 2 d ~ 2 . 

B3: Given the bounds on vol(Y) and vol(AT) < 2 d_1 , Lemma [ 4 ] tells us that the edges of X have 
axes from |cr x (l)|,..., \a x (d — 1)|, and the edges of Y have axes from |tr y (l)|,..., \<J y (d — 3)|. 
Now let a = \a x (d)\. By Theorem |TT| depth(cr y ,a) < depth((7 x , a) + 1 = 1 and therefore a is 
not included in |<r w (l)|,..., \(T y {d — 3)|. If a = c, it follows that X and Y lie in half-cubes 
that together constitute a unit cube of volume 2 d , as in case Bl. 

Otherwise, if a c, it follows that E may contain multiple edges of direction c but does 
not include any edge with direction a. Therefore E lies completely in a box that spans 
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two 1-cubes in dimension c, half a 1-cube in dimension a, and one 1-cube in the remaining 
dimensions. The volume of this box is 2 d . 


C: By Lemma 12 each set of vol(-E) — 1 edges of Ak is contained in a unit cube of |"log 2 (vol(£ ; ) — 
1)] + 1 = Llog 2 (vol(.E) — 2)J + 2 dimensions, of volume at most 4(vol(.E) — 2). 

□ 


4 General construction method in three and more dimen¬ 
sions 


4.1 Extended curves and local edge distance 

In Section [2] Theorem [7J we learned about suf¬ 
ficient and necessary conditions for well-folded 
curves in general, and in Section [3j Theorem |11[ 
we learned about specific conditions for hyper- 
orthogonal well-folded curves. It remains to 
show that curves satisfying both the general 
and the specific conditions actually exist. 

In this section we will combine the condi¬ 
tions of Theorems [7] and Qj] to derive conditions 
on the entry and exit points and the isometries 
used in the construction of hyperorthogonal 
well-folded curves. We will show how to con¬ 
struct curves that satisfy all conditions, for any 
d > 3 (recall that for d = 2, we have Hilbert’s 
curve). 
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Figure 3: in black: G( 3) with the directions of 
its edges; in grey: an extension of G(3) with an 
entry edge (d) and an exit edge (— (d— 1)), with 
the edge distance table for each vertex according 
to Definition fTU 


Definition 14. The edge distance of the axis 
a £ {1,..., d} to the vertex v within the curve 
C, denoted ed(C, v, a), is the distance along C 
between v and the closest edge with axis a; more 
precisely, ed (C,v,a) is one less than the length 
of the smallest subcurve of C that includes v 
and an edge with axis a. (For a small example, 
see Figure [3p 

Theorem [TT| has a remarkable consequence: 

Lemma 15. Let Ao ,... A^+i be a sequence of well-folded hyperorthogonal curves constructed by 
inflation. Then we have, for all axes a £ {l,...,d} and all vertices W; of Ak, depth(<Ji,a) < 
ed (A k ,Vi,a). 

Proof. The proof goes by induction on increasing edge distance. If ed (Ak,Vi,a) = 0, then a € 
{|ej|, |e,_i|}, and, by the first condition of Theorem 11 we have depth(cr, ; , a) = 0 = ed (Ak,Vi,a). 


Now suppose ed (Ak, u,, a) > 0. Then we can choose j £ {i — 1, i + 1} such that ed(4fc, Vj, a) = 
ed(4/ c , Vi, a) — 1, and by induction we can assume depth(c7j, a) < ed (Ak,Vj,a). Then it follows from 


the second condition of Theorem 11 that we have depth(crj, a) < depth(oj, a) + 1 < ed (Ak,Vj,a) + 
1 = ed(A k ,Vi,a). □ 


Lemma 15 gives us the following idea for an algorithm to specify the permutations a., |, except 
for the order of the last two elements: simply sort all axes a £ {1,..., d} by order of decreasing 
edge distance ed (Ak,Vi,a). In fact, as we will show, a version of this algorithm suffices, that only 
considers edge distances within small subcurves. For this purpose we define the notion of extended 
curves, which can be seen as curves together with an indication of how the curve is connected to 
preceding an succeeding curves: 
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Definition 16. An extended curve is a curve that is extended with an entry edge leading to the 
first vertex (the entry point) and an exit edge originating from the last vertex (the exit point). The 
origin of the entry edge and the destination of the exit edge are not considered to he part of the 
curve. 


We use prime symbols to distinguish extended curves from non-extended curves: when B is a 
curve, we may use B' to denote a particular extension of B, and when B' is an extended curve, we 
use B te denote the curve without the extensions. Note that by our definition, B and B' always 
have the same vertices; they only differ in the number of edges. In particular, if we extend an 
approximating curve A k that has edges (e^i,... ,ek,K- 1 ), we denote the extended curve by A ' k , 
the entry edge by e^.o and the exit edge by ek t K■ If B is a subcurve of A', then the entry edge of 
B' is the edge that leads to the entry vertex of B in A', and the exit edge of B' is the edge that 
originates from the exit vertex of B in A'. In particular, the extended child curve C k r of A' k would 
be Ckj. extended with entry edge (ek,i-i) and exit edge ( ek,i)■ 

The definition of well-foldedness (Definition [5]) can be applied to extended curves, with the 
base case that an extended curve that consists of only an entry edge, a single vertex, and an exit 
edge, is well-folded. The conditions for well-foldedness from Theorem [7] are applicable as well. In 
that case it is natural to require that we would obtain a valid curve if we would add the origin of 
the entry edge and the destination of the exit edge to the curve. This can be ensured as follows: 
we take the entry edge into account by extending the second condition of Theorem [7] to the case 
i = 0; given the first condition, the second condition can also be written as: sign (er” 1 ^)) = 1 
if and only if \ef = \a t (d)\, and we take the exit edge into account by extending this form of the 
condition to the case i = K. 

The definition and conditions of hyperorthogonality (Definition [9] and Theorem 11) can be 
applied to extended curves, if, in condition 1 of Theorem m we also take the entry and exit 
edge into account. Concretely, this means condition 1 should be extended with depth(cr 1 , eo) = 0 
and depth((7if, ex) = 0. The definition of edge distance (Definition 14) and its relation to 
hyperorthogonality (Lemma 15) can now be applied directly to extended curves. 

We can now define a version of edge distance that only considers small subcurves: 


Definition 17. Let A' k he an extended well-folded curve obtained by inflation from A' k _ l . Let v he 
a vertex of Ak, let a £ {1,... ,d} be any axis, let the subcurve C of Ak be the child curve of A , k _ 1 
that contains v, and let C be the extension of C within A' k . We define the local edge distance of 
the axis a to the vertex v within the curve A' k , denoted led(A^., v, a), as ed (C',v,a). 


4.2 Hyperorthogonal curves from inflation of extended curves 


Lemma 18. Suppose we construct a sequence of extended well-folded curves A ' 0 , A[,... by inflation 
such that the elements of each permutation \<Jk,i\ are sorted by order of decreasing local edge distance 
to Vk,i in A k . Then these permutations satisfy conditions 1 and 2 of Theorem 11 


Proof. The proof goes by induction on increasing k. As the base case we take k = 0, and observe 
that Aq, which contains only a single vertex and two edges, trivially satisfies Theorem 11 Now 


suppose k > 0 and the permutations associated with the vertices of A' k _ l satisfy the conditions of 
Theorem We will now show that, if we choose the permutations 01 , 02 , ■■ ■ associated with the 
vertices v\, v%,... of A' k in such a way that the elements of \oi\ are sorted by order of decreasing 
local edge distance to Vi in A), then these permutations satisfy the conditions of Theorem 11 as 
well. 


Consider any vertex Vi in A' k . 


Let C' be the extended child curve of A' k _ 1 that contains vu 


and let r be the signed permutation such that C is a translate of r(G(d). Since C' includes both 

edges of A' k that are incident on 1 we have led(A' fc , Vi, a) = 0 if and only if a £ {|e*_i|, |e*|}. 

Hence, these two axes will be placed at the last positions within |cq|, so that depth(cri,a) = 0, and 
condition 1 of Theorem [lT] is satisfied. 

For condition 2, observe that in C, being a transformation of G(d ), the edges with axis |r(l)| 
and edges with other axes alternate, starting and ending with an edge with axis |r(l)|. By the 
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induction hypothesis, A' k _ 1 satisfies the conditions of Theorem I 111 which implies that the edges 
immediately preceding and following C in A' k have axes with depth zero in r. Therefore these 
axes differ from |t( 1 )|, which has depth d — 2 (remember that this section is concerned with 
d-dimensional curves for d > 3). Thus, also in C' edges with axis |r(l)| and edges with other axes 
alternate. 

Now suppose, for the sake of contradiction, that there are two axes a ^ b, both different from 
|r(l)|, such that led(H' fc , v t , a) = led(H(,, v t , b) for some Vi £ C'. Then C' must contain an edge 
sequence of even length, more precisely 2 * led(H' fc , Uj, a) + 2 , with Vi in the middle, that starts with 
an edge with axis a and ends with an edge with axis b. However, this contradicts the fact that 
edges with axis |t( 1)| and edges with other axes alternate. Hence, apart from a pair of axes with 
local edge distance zero (among which |r(l)|), no pair of axes has the same local edge distance. 

If we increase i by one while staying inside the same child curve C, each local edge distance 
changes by at most one, and therefore each axis can move up or down in the order of r by at most 
one position. This establishes condition 2 of Theorem [TI] as long as we stay in the same child curve, 
that is, as long as \i/2 d ] does not change, that is, for all* ^ 0 (mod 2 d ). 

If i = 0 (mod 2 d ), we need to take more care, as Vi lies in a child curve Cj of At* while Wj+i 
lies in the next child curve Cj+ 1 . Now, since A' k _ 1 satisfies condition 1 of Theorem fTTl there must 
be g, h £ {d — 1 ,d} such that |cr fc _i ^ (< 7 ) | = \crk-i,j+i{h)\ = \ek~ij\ where e^-ij isthe direction 
ek t % of the edge (uj, Uj + i) in A' k . Define g', h! £ {d — 1 , d} by g' 7 ^ g and h! ^ h. 

Now, for Vi, sorting by decreasing edge distance within 6 " results in |cr fc)i | = 

[Wk-i,j(9% Wk-i,j(d - 2)|, |- 3)|,..., |crfc— 1 (2)|,K-ij(l)|, \a k - 1}i {g)\] , 

where the order of the last two elements is undetermined, and likewise 

for Vi+ 1 , sorting by decreasing edge distance within C' +1 results in \(Jk,i+i \ = 

(d ) |, | {d 2) |, (d 3)|,..., \(Jk—\ 7 j-\-i (2)|, | (Jk-i 7 j-\-i (1) |, lu'/c— 1 , j+iWI] i 

where also the order of the last two elements is undetermined. For a = \ek-i,j\ we thus have 
depth(o’fc j i(a)) = depth(crfc i j_ ) _i(a)) = 0 , and for a ^ \ek-i,j\ we have depth(crfc : j(o)) = d — 2 — 
depth(<Tfc_i j;; -(a)) and depth(cr/ Ci i + i(a)) = d — 2 — depth(crfc_i J+ i(a)). By the induction hypothesis, 
and CTfc-ij+i satisfy condition 2 of Theorem [IT] and therefore we have | depth(crfc_i j (o)) — 
depth(crfc_i J+ i(a))| < 1 , and hence | depth(crfc i j(a)) — depth(crfc,i+i(a))| < 1 , which establishes 
condition 2 of Theorem [TlJ O 

The above lemma still leaves the order of the last two elements of each |< 7 j| undetermined, since 
these are always the two axes with edge distance zero. To prove that hyperorthogonal well-folded 
curves exist, it now suffices to show that we can order the last two elements and choose the signs 
of each ai such that the conditions of Theorem [7] are satisfied. We obtain: 

Theorem 19. For each choice of an entry direction eo.o and an exit direction eop and for each 
choice for the signs of a k \(j) for all k and j such that sign(cr^T^(efc j o)) = 1 for all k, there is a 
unique hyperorthogonal , well-folded space-filling curve f approximated by A' 0 , A[,..where each 
curve A! k with k > 0 is constructed by inflation from A! k _ x and the elements of each permutation 
|(j k t i | are sorted by order of decreasing local edge distance to Vi in A' k . 

Proof. For each level k, we generate A' k as follows. We loop over all i £ {1,..., K — 1}, where 
K = 2 d * k , and proceed as follows. The conditions of Theorem m require sign(cr“ f _ 1 1 (ei)) = 1. We 
now choose |ui(d)| such that Icr^d)! = |e*| if and only if sign(cr“ T (e i )) = 1 : this is always possible 
since || is among the last two elements of |(X;| whose order was undetermined. Thus we satisfy 
the second condition of Theorem [7] for j = \d\. With |cr,| completely determined, we can now fill 
in the remaining signs of 1 such that they fulfill the first condition of Theorem [7J Finally, we 
determine |cr/c(c?)| as dictated by the exit direction ej<- in the same way as we determined |cq(ci)| 
for i < K. □ 

Note that if A' 0 , A\,... is a set of extended hyperorthogonal well-folded curves constructed by 
inflation, then they are approximating curves of a hyperorthogonal well-folded space-filling curve. 
Note, however, that not every hyperorthogonal well-folded space-filling curve can be described by 
such a set A’ {) , A\,.... There can also be hyperorthogonal well-folded space-filling curves that can 
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Figure 4: Each cube in this figure shows a grey curve serving as a very rough sketch of Ak, with 
its entry and exit on the interior of two non-opposite (d — 1)-dimensional faces of the cube, (a) It 
is easy to connect up four copies of such a curve, (b) Assembling more than four copies requires 
rotating the cubes to bend the path into other directions. In general, the rotations will break the 
connections between one copy and the next. 


be described by a set of non-extended hyperorthogonal well-folded approximating curves Aq,Ai,... 
which cannot be extended to a set A' 0 , A ' 1 ,... of hyperorthogonal curves in such a way that each 
curve A! k is equal to Ak extended with an entry edge (eo,o) and an exit edge (eo,i). Examples 
would include symmetric space-filling curves, closed space-filling curves (that is, curves that start 
end end in the same point), and space-filling curves that start in the interior of the unit cube. 


5 Self-similar curves in three and more dimensions 


5.1 The challenge 

By Observation |8j a choice of signs of (j) for all k and j specifies the starting point /(0) of 
the space-filling curve / in Theorem |19| Thus, the proof of Theorem 19 is a constructive proof 
that a hyperorthogonal, well-folded space-filling curve exists for any choice of starting point on the 
boundary of the unit hypercube. 

In a practical setting, such as described in Section |1.1[ one may want to sort points in the order 
in which they appear along the curve. To this end we need a comparison operator that decides 
which of any two given points p and q comes first along the curve. We can do so by determining the 
largest k such that there is a hypercube Hk,i, corresponding to a vertex Vk,i of Ak, which contains 
both points. Then we can use cr^j to determine in which order the 2 d subcubes of this hypercube 
are traversed, and in particular, in which order this traversal visits the two subcubes containing p 
and q. The efficiency of the comparison operator now depends on how efficiently we can determine 
CTk.i for any k and i. Unfortunately, straightforward application of Theorem would require us to 
derive in an incremental fashion that explicitly constructs all ak,j for all j < i. In practice 
we will need a less time-consuming way to derive <Jk,i- This seems rather difficult to achieve if we 
allow ourselves to choose the signs of the permutations <JkA arbitrarily. 

To enable us to determine a permutation more efficiently, we will, in this section, restrict 
the curves to be self-similar , that is, any approximating curve Ak +1 is the concatenation of 2 d 
isometric copies of Ak- Recall that taking the reverse is also an isometric mapping. 

Note that for d = 2, Hilbert’s original curve is the only self-similar well-folded curve. So for the 
purposes of Sections [5 . 2| to [576] we will assume d > 3. 

In Section [572] we find that all self-similar hyperorthogonal well-folded space-filling curves fit 


the framework of Theorem 19 that is, they can be described by a series of extended approximating 
curves. Moreover, we find that the only extended 1-curves that are relevant for the study of self¬ 
similar hyperorthogonal well-folded space-filling curves are isometries of one particular extension 


of G{d). As shown in Section 4.2, a (not necesssarily self-similar) hyperorthogonal well-folded 
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space-filling curve does not need to start at a corner of the unit cube, but may start on the interior 
of a (d — l)-dimensional face. In Section 5.3 we set up notation which is helpful in distinguishing 


different possible locations for entries and exits. In Section pci] we analyse how the choice of the 
entry of Cpi propagates to the other child curves C \$,..., Ci t n of A±. We will find that if the 
entry point of A± lies on the interior of a (d — l)-dimensional face, the exit point will also lie on 
the interior of a face. If the entry and exit of an approximating curve Ak indeed lie on different 
but non-opposite faces, it is now trivial to connect up four copies of A k to make a cycle using only 
reflection and reversal transformations, see Figure [4^ a). 

However, to get beyond four copies of Ak and assemble 2 d copies to make a full Ak+ 1 , we need 
to rotate some copies of Ak in various ways to bend the path into all d dimensions. The difficulty 
is to ensure that despite the various rotations, the connection points on the faces will still match 
up. This will not automatically be the case, see Figure [4^b) . As we will find in Section 5.5 this 
forces most of the coordinates of the entry and exit points to be equal, so that these points become 
invariant under the necessary transformations. 

This strongly restricts the possible shapes of self-similar hyperorthogonal well-folded curves, 
but not too much: in Section 5.6 we find that such curves do in fact exist. It turns out that for 


any d > 3, only two different starting points (modulo rotation and reflection) exist for such curves. 


5.2 Extensions in self-similar curves 

As noted in Section |4.2[ not all hyperorthogonal, well-folded space-filling curves may be approxi¬ 
mated by extended hyperorthogonal, well-folded curves. However, for self-similar curves this can 
always be done. But before proving this, we will first have a look at extensions of G(d). In extended 
hyperorthogonal well-folded approximating curves, we find only one particular extension of G(d): 

Definition 20. Let G'{d) be the concatenation of (d), G(d), and (—(d — 1)). 

Lemma 21. Let A ' 0 ,..., A' k , 1 be a sequence of extended hyperorthogonal well-folded curves con¬ 
structed by inflation. Then each extended child curve C' ki is the image of an isometry ofG'(d). 


the axes of e,_i and ej in C[ must be |afld — 1)| and \ai(d)\. This leaves two possibilities for 
matching axes to edges. The first possibility is to put the edge with axis cr,;(fi) at the beginning 
and the edge with axis |er;(ci — 1)| at the end; the signs of the directions follow from the fact that 
the vertex preceding Ci in Ak+\ and the vertex following Ci in Ak+\ must lie outside Ci. The 
second possibility is to put the edge with axis |crfld — 1)| at the beginning and the edge with axis 
|cr,;(d)| at the end, so that we obtain a concatenation of {afld — 1)), o' i (G(d)) and {afld)), which is 
the reverse of aflG'{d )) with reflection in coordinate ofld). □ 

Note that the proof of Lemma [21] does not require the approximated space-filling curve to be 
self-similar. 

Definition 22. If, in the above lemma , the isometry that maps G' (d) to C[ is composed exclusively 
of rotation, reflection, and translation, then we say C' k i is of type 0; otherwise, that is, if the 
isometry that maps G'(d) to C' k t involves reversing the curve, then we say C k i is of type 1. We 
denote the type of C’ ki by Tk 

As always, when the first subscript to T is clear from the context, or when a statement holds 
for any value of the first subscript, we may omit the subscript. 

In the following observation we use the Iverson bracket notation: when P is a expression that 
evaluates to true or false, then [P] = 0 if P is false, and [P] = 1 if P is true. The observation is 
the following: the type of an isometric image C[ of G'(d) is zero if and only if |cq(d)| is the axis 
of the entry edge e»_i, and the type is one if and only if |cq(d)| is the axis of the exit edge e,. In 
other words: 

Observation 23. Ti = [|cq(d)| |e*_i|] = [|cr i (d)| = |e,|]. 


Proof. Consider an extended child curve C[ of A 


k- 


Since A' k is hyperorthogonal, by Theorem 
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We can now prove that all self-similar hyperorthogonal, well-folded space-filling curves can be 
approximated by a series of extended hyperorthogonal, well-folded curves: 

Lemma 24. Let f be a self-similar, hyperorthogonal , well-folded space-filling curve. Then an 
isometric copy of f is approximated by a series of extended hyperorthogonal, well-folded curves 
Aq, A[,..., where Aq = {d,—(d— 1)), <To,i = [l,...,d], and each extended curve A' k with k > 0 is 
obtained by inflation from A' k _ 1 . 

Proof. For any k > 0 and 1 < i < 2 d , let B\ i^ be the fc-curve that is a subcurve of the non- 
extended approximating curve A k+ i and results from k steps of inflation of the vertex Vi of A\. 
Since A^+i is hyperorthogonal and well-folded, the extended curve B[ t k that consists of B\^^ 
with entry edge (ej_i) and exit edge {ef) must also be hyperorthogonal and well-folded. It follows 
that B[ i0 , B[ 11 , B[ i2 , ■ ■ ■ is a sequence of extended hyperorthogonal, well-folded curves that 
approximate the space-filling curve fl which consists of / restricted to the hypercube Hi that 
corresponds to u t . By Lemma |21[ it follows that there is an isometric transformation that maps 
B[ i l to G'(d), and thus, B' li0 to (d,—(d— 1)). 

Because / is self-similar, the same series of curves that approximates /,; also approximates /, 
up to isometric transformations. □ 

Lemma 25. Let F 1 be an extended hyperorthogonal well-folded curve obtained by one step of 
inflation from G'(d), and let R! be an extended hyperorthogonal well-folded curve obtained by one 
step of inflation from G'{d). Then no non-reverse isometry of R! can visit its vertices in the same 
order as F'. 

Proof. Suppose r is a non-reverse isometry, expressed by a signed permutation, such that t(R') 
visits its vertices in the same order as F', and hence the axes of the edges of t(R') and F' are the 
same, apart from, possibly the entry and the exit edge. Then r(G 7 ((i)) must also visit its vertices 
in the same order as G'(d), so r = [1,..., d — 1, —d], and |r| is the identity permutation. However, 
since the entry edge of t(R') has axis d — 1 while the entry edge of F' has axis d, the axes of the 

edges of t(R') and F' must differ in the first child curve of r(G"(d)) and G'(d), respectively. □ 

Note that another way to put the last line of the lemma is to say that any non-reverse isometry 
of F' must differ from R' in more than just the entry and/or exit edge. 

Corollary 26. For d > 3, any d-dimensional self-similar hyperorthogonal well-folded space-filling 
curve is asymmetric. 

5.3 Relative coordinates of entries and exits 

In the following subsections, the following notation will be helpful. 

Definition 27. Let entfc im , extfc jm : {l,...,d} —> {0,..., 2 k+1 — 1} be functions that give the 
coordinates of the entry and exit point of Ck, m , that is, the entry point of Ck, m has coordinates 
( ent m (l),..., ent m (d)) and the exit point has coordinates (ext m (l),..., ext m (d)). 

Note that C/ im is a 1-curve that is a subcurve of A k+ i, we have ent k,m(j) = isneg^/'^/j)) 
(mod 2), and ext k,m(j) = ent k,m(j) (mod 2) if and only if |<x“j n (j)| ± d. 

Definition 28. The relative coordinate vector of a vertex v is the vector r such that r[j] = 0 if 
v\j\ mod 4 £ {0, 3}, and r\j] = 1 if v[j] mod 4 £ {1, 2}. 

The relative coordinates of a vertex v n of Ak+\ tell us, for each dimension, whether the vertex 
is on the outside (0) or on the inside (1) with respect to the 2-cube that results from inflating the 
inflation of the vertex Vj of A^- 1 , where j = \n/D 2 ~\. 

Definition 29. Let rlentfc )m ,rlextfc iTO : {1,... ,d} —¥ {0,1} be functions that give us the relative 
coordinates of the entry and exit point of Ck,m- 
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Observation 30. 

rlent fc>m (j) = (ent k , m (j) + v k ,m[j\) mod 2 , or equivalently, 
rlent k,m{j) = (isneg(<T^(j)) +v k , m \j}) mod 2 , and 
rlext & ,m (j) = (extfc ,ra GO + vk ,m [j]) mod 2. 

Note that in the above observation, entfc jm (j) is a coordinate of the entry of Ck,rm which is a 
vertex of Ak+ 1 , while Vk, m [j] is a coordinate of a vertex of Ak . In fact, Vk, m [j] = L en tfc,mG)/2j 
and entfc jm (j) must be either 2 * Vk, m [j] or 2 * Vk, m [j\ + 1- Thus, if rlentfc im and Vk, m are given, 
this determines entfc jm and hence, the signs of cr^T m . 

Observation 31. 

rlext k,m{j) = rlent fejm (j) + [\cr k ,m(d)\ = j] (mod 2); 
for m < D k we have rlextfe jm = rlentfc im+ i. 

As always, when the first subscript to ent, ext, rlent or rlext is clear from the context, or when 
a statement holds for any value of the first subscript, we may omit the subscript. 


5.4 Relation between entry and exit of a 2-curve 

A direct consequence of Lemma |24| is that for a self-similar curve we may assume, without loss 
of generality (modulo reflection, rotation and reversal), that A! x = C 01 = G'(d ), so with type 
Top = 0, entry edge (d) and exit edge (— (d — 1)). Moreover, in A 2 , we should have ^i[d] = 0 
and vx[d — 1] = 0, where K = D 2 = (2 d ) 2 , so that the child curves Cpq and . o m A 2 can be 
extended with, respectively, the same entry edge (d) and the same exit edge (—(d — 1)) as A±. Note 
that we can rewrite the conditions on vi[d] and vx[d — 1] as rlentiq(d) = 0 and rlexti j £)(d— 1) = 0. 

In this subsection we consider extended hyperorthogonal well-folded approximating curves A\ 
and A 2 that fulfill these basic conditions, that is, A\ = G'(d), rlentiq(d) = 0 and rlexti ,D(d) = 0, 
without assuming, at this point, that A\ and A'. 2 are indeed approximations of a self-similar 
space-filling curve. In particular, in this subsection we analyse how the choice of the entry point 
of Cyi propagates to the other child curves C\ i2 ,... ,C\^d of A\. Because the whole subsection 
focuses on the child curves of A[ that constitute A ' 2 , we will omit the first subscripts to C , a , T, 
rlent and rlext: they would always be 1. 

Definition 32. Let ui be the permutation [d — 1, 2, ..., d — 2, d, 1] . 


By tracing the relative coordinates of the entry and exit points through the child curves of 
A\ that make up A 2 , using the conditions of Theorems [7] and 11 we find rlext£> = rlenti o ui (in 
Lemma 36). To prove this we need the following three lemmas from which the proofs can be 
skipped at first reading. 


Lemma 33. 


• rlenti (d) = 0. 

• If rlenti(l) = 0, then |<Ji(d— 1)| = d and |cri(d)| = 1, 
otherwise |<7i(d— 1)| = 1 and |cri(d)| = d. 

• Ti = 1 — rlenti (1). 


Proof. The first item, rlenti(d) = 0, follows from the fact that the entry edge is (d), by the 
assumptions of this subsection. The second item follows from the fact that we need rlexti (1) = 1 to 
be able to connect C\ to C 2 with the first edge of G(d), which has direction 1. The third item follows 
23 by Ti = [|<Ji(d)| = |ei|] = [|oi(d)| = l] = [rlenti(l) = 0] = l-rlenti(l). □ 


from Observation 


Lemma 34. For even i <2 d 


have: 


• rlenti(l) = 1; for 1 < j < d, rlenti(j) = rlenti(j); rlentj(d) = rlenti(l). 
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• //rlent 2 (|ei|) =0 then |<7j(d— 1 )| = |cr;_|_i(d— 1 )| = 1 and |<7j(d)| = |<Ji+i(d)| = \ei\, otherwise 
\a t (d- 1 )| = |cr i+ i(rf — 1 )| = |e*| and |<+(d)| = |cr i+ i(d)| = 1 . 

• T i+ i = rlent 2 (|ei|); T) = 1 - T i+1 . 


Proof. Throughout this proof, all additions and subtractions are to be interpreted modulo 2. 

We first handle the case i = 2. 

It is straightforward to calculate the relative entry function rlent2 from rlenti and <7i (d) 
using Observation |31| In particular, with the second item of Lemma |33| we get rlent2(l) = 
rlenti(l) + [|oi(d)| = l] = rlenti(1) + (1 — rlentr(l)) = 1; for 1 < j < d we have rlent 2 (j) = 
rlenti (j) + [|<ri(d)| = j\ = rlenti (j); and rlent 2 (d) = rlenti (d) + [|cri(d)| = d\ = 0 + rlenti(l). 

Now, because i is even, by Lemma|2]we have |e*| ^ 1 and |ei_i| = |e i+ i| = 1. With respect to 
axis Cj|, the exit point of C* must be on the inside of the 2-cube traversed by A' 2 , otherwise it 
cannot be connected to the next child curve Ci+± by an edg e with axis |e» |. In other words: we 
must have rlexti(|ej|) = 1, and therefore, by Observation 31 [|cr,;(d)| = |e*|] = 1 — rlenti(|ej|). 
Therefore, if rlentj(|e,|) = 0, then |cr,(d)| = | eff and therefore, by Theorem 

| &i{d — 1)| = |e*_i| = 1; moreover, rlentj + i(l) = rlentj(l) + [|<7i(d)| = l] = rlenti(1) = 1, 
[|<Tj+i(d)| = l] = 1 — rlentj+i(l) = 0, in other words, |<7i+i(d)| ^ 1 = |ej+i|, and, by Theorem 
\<j i+1 (d)\ = |ei| and |cr i+1 (d- 1)| = 1. 

Otherwise, if rlenti(|ej|) = 1, then \(Ji{d)\ |ei|, so, by Theorem 11 |crj(d)| = 1 and |aj(d—1)| = |ej|; 

moreover, rlent, + i(l) = rlenti(l) + [|tTi(d)| = l] = rlenti(l) + 1 = 0, so [|<T, + i(d)| = l] = 
1 — rlenti + i(l) = 1, in other words, |(7i + i(d)| = 1 = |ei+i|, and, by Theorem [II] |cri + i(d— 1)| = \d\. 

Ti +i = [ki+i(d)| ^ |ei|] = 


11 


so 


11 
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The third item of the lemma now follows from Observation 
rlenti(|ei|), and Tj = [|<Ti(d)| = |ei|] = 1 - rlenti (|ei|) = 1 - T i+1 . 

The cases i > 2, for even i, then follow by induction, using that <Ti + i(d) = <Ji(d), so that, by 
Observation 31 rlenti+2 = rlext»+i = rlenti. D 


Lemma 35. 


• rlentu(l) = 1; for 1 < j < d, rlentu(j) = rlenti(j); rlent£>(d) = rlenti(l). 

• //rlentu(d — 1) = 0 then \(To(d — 1)| = d — 1 and <Tu(d) = —1, 
otherwise \ao{d — 1)| = 1 and ao{d) = —(d — 1). 

• Tq = rlentu(d — 1) = rlenti (d — 1). 

Proof. The first item is actually proven in the last line of the proof of Lemma [34] with * = D — 2. 
The second and third item follow from straightforward calculations, similar to those of the previous 
lemmas, where we use the fact that, by the assumptions of this subsection, we have eo-i — (—1) 
and eo = (—(d— 1)), and therefore rlentu(l) = 1 and rlext_o(d — 1) = 0. □ 

Lemma 36. rlext /) = rlenti o uj. 

Proof. Straightforward rewriting of the equations in Lemma [35] using Observation |31[ yields: 

• rlextu(l) = rlentu(l) + [|cr_D(d)| = l] mod 2 = rlentu(l) + [rlentu(d — 1) = 0] mod 2 = 

rlentu(l) + [rlentu(d — 1) = 0] mod 2 = 1 + (rlenti (d — 1) + 1) mod 2 = rlenti(d — 1); 

• for 1 < * < d — 1, we have rlext£>(?) = rlenti (*) + [|cr£) (d) | = *] mod 2 = rlenti (i) + 0 = 

rlenti (i); 

• by the assumptions of this subsection, rlextu(d — 1) = 0, which equals rlenti (d); 

• rlextu(d) = rlentu(d) + [|cr£>(d)| = i] mod 2 = rlenti(l) + 0 = rlenti(l). 

This establishes rlextu = rlenti ou> with w as in Definition [32] □ 


20 









5.5 Possible entry points of self-similar curves 


In this subsection we will first use the similarity between the 2-curves that make up A 3 to prove 
Lemma 42 which states that rlenti(j) should be the same for all j £ {1 ,d — 1}. After that, we 
will use the similarity between A 2 and the 2-curve that forms the beginning of any approximating 
curve Ak ( k > 2), to prove Lemma 43 which states that rlentfc = rlenti for all k > 1. From that we 
will derive Theorem |45[ which essentially says that for any fixed d, there are only two points that 
may be the starting point of a d-dimensional self-similar, hyperorthogonal, well-folded space-filling 
curve. 

Let A! x , A' 2 , A '3 be extended hyperorthogonal well-folded approximating curves of a self-similar 


space-filling curve, fulfilling the assumptions which we made, without loss of generality, in Section 5.4 
When we inflate A ' 2 to obtain A' 3 , so that a 2-curve replaces each vertex of A\. the relative 
coordinates of each 2-curve’s exit point should equal the relative coordinates of the next 2-curve’s 
entry point—otherwise the 2-curves would not be connected by an edge. 


Observation 37. Because of self-similarity, the 2-curve replacing Vi of Ai must itself be an 
non-reverse isometry of either A 2 if Ti = 0, or A 2 ifTi = 1. 

Note the either-or in the above observation: by Lemma |25| the 2-curve replacing Vi cannot be 
a non-reverse isometry of both A -2 and A 2 at the same time. 

As a result of the transformation cy*_ 1 , the relative coordinates of the exit point of the 2-curve 

replacing Vi -1 of A\ are given by the function rlenti ocuo of]_ if Tj_ 1 = 0, and by rlenti o \<j~\\ 
if 1 = 1. The relative coordinates of the entry point of the 2-curve replacing Vi are given by the 
function rlenti o |cr^ 1 1 if Tj = 0, and by rlenti oojo |cr~ 1 1 if Tj = 1. Thus we get: 

Lemma 38. 


• IfTi -1 = 0 and Ti = 0, we have rlenti owo \cr~\\ = rlenti 0 \&f 1 \ 

• IfTi -1 = 0 and Ti = 1, we have rlenti owo Icr^l-d = rlenti owo \crf 1 \ 

• IfTi -1 = 1 and Ti = 0 , we have rlenti o \a~\\ = rlenti 0 Wff \ 

• IfTi_i = 1 and Ti = 1, we have rlenti o \a~} 1 \ = rlenti °wo \vf 1 \ 
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which states 


We will now analyse the possible successions of types Ti and permutations &i for the vertices Vi of 
Ai, where i £ {1,..., 2 d }. We will do so in four lemmas, concluding with Lemma 
that rlenti (j) should be the same for all j £ {1,..., d — 1}. 

Lemma 39. There is a j £ { 2, 3,..., D} such that Tj = Tj_\. 

Proof. There is an even i < D = 2 d (specifically, we may choose i = D/A or i = 3D/4) such that 
\ei\ = d — 1. By Lemma 34 we have T i+ 1 = rlenti(d — 1), which, by Lemma 35 equals Tj-j . Thus 


the sequence Tj + 1 , Tj_|_ 2 , ..., Tp consists of an even number of types where the last equals the first. 
This implies that a strictly alternating type sequence is not possible. □ 


Lemma 40. rlenti (1) = rlenti (d — 1). 


Proof. Let i be the largest i from {2,3,..., D} such that Tj = Tj_i (there is always such an i, by 
Lemma [39| . We distinguish three cases: (i) |crj_i(l)| = |<jj(l)|; (ii) |crj_i(l)| ^ |cq(l)| and d > 4; 
(iii) l<A_i(l)l M 1 )! and d = 3. 


In the first case, let x be crj(l), so we have \a~\{x)\ = |cr“ 1 (a;)| = 1. Then, from Lemma 38 
evaluating the functions on both sides for x, we find, both in the case of Tj = Tj_i = 0 and the 
case of Tj = Tj_i = 1, the following: rlenti(l) = rlenti(w(l)) = rlenti(ci — 1). 


In the second case, we have, by Theorem 11 |cr. i _i(1)| = |<7j(2)| and |(jj_i(2)| = |oj(l)|. Let 


x = \<Ti(l)\;y = |crj(2)|, so we have |cr i _ 1 1 (y)| = \a i 1 (a;)| = 1 and |cr i _ 1 (x)| = \a. l 1 (j/)| = 2. Since 


d > 4 we have ui(2) = 2. Now, if Tj_i = Tj = 0, Lemma 38 gives us rlenti(w(|c7j_ 1 (a;)|)) = 


lenti(|<Tj i (®)|) <-> rlenti(2) = rlenti(1) and rlenti(w(|CTj_ 1 (y)|)) = rlentidcq (y )|) O rlenti(d — 
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1) = rlenti (2), so rlenti(l) = rlenti(2) = rlenti(d — 1). Otherwise we must have Tj_i = Tj = 1 


and Lemma 38 gives us rlent 1 (|tr-_ 1 (x)|) = rlenti(w(|cr i i (a;)|)) O rlenti(2) = rlenti (d — 1) and 
rlent i(|<t“ 1 1 (i/)|) = rlenti(w(|cr“ (y)|)) o- rlenti(l) = rlenti( 2 ), so, again, rlenti(l) = rlenti( 2 ) = 
rlenti (d — 1 ). 

The third case does not occur, since for d = 3, the proof of Lemma [39| yields T s = TV, so i = 8 . 
Since G'(3) ends with (..., —2, —1, —2), both Vj and Vg are incident on edges with axes 1 and 2, 
and both <77 and ag must have the remaining axis, 3, at depth 1, thus |< 7 j_i(l)| = |< 7 j(l)| = 3. O 

Lemma 41. rlenti (j) = rlenti (j — 1) for all j £ {3,4,..., d — 2}. 

Proof. By Theorem [TT] we must have depth(<7/i, a;) = 0 for some h and x = |<7i(l)|, so depth(<Ji, x) = 
d— 2. Since depth differs by at most one between successive permutations < 7 j, there must be, for any 
j £ {2, ... ,d — 2}, an i such that \cn(j)\ = — 1)| = x. Note also that for j £ {2,..., d — 2}, 

we have u>(j ) = j. Hence, from Lemma 38 evaluating the functions on both sides for x, we find 
rlenti {j - 1) = rlenti (j). □ 

Lemma 42. rlenti (j) = rlenti (J — 1) for all j £ {2,..., d — 1}. 

Proof. If d = 3, the lemma is equivalent to Lemma [40j Otherwise, choose i and x such that 
|< 7 j( 2 )| = |cq_i(l)| = x (such i and x exist, as observed in the proof of Lemma[4l|). Now, if Tj = Tj_i, 
the proof of the second case of Lemma 40 tells us that rlenti (1) = rlenti(2) = rlenti (d — 1), and 
Lemma |42| follows by combining this fact with Lemma |4lJ It remains to discuss the cases in which 
T t ^T z _ 1. ' _ 

we have rlenti(d— 1 ) = rlenti(w(l)) = rlenti(o;(cr i _\(a;))) = 

and 
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If Tj_i = 0 and Tj = 1, by Lemma 

rlenti(^(cr .” 1 (a;))) = rlenti(w(2)) = rlenti(2). Combining this with Lemmas 


rlenti(j) = rlenti(j — 1) for all j £ {2,3,..., d — 1}. 


40 


41 


establishes 


If Tj_i = 1 and Tj = 0, by Lemma 
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we have rlenti ( 1 ) = rlenti(cr i _ 1 (x)) = rlenti(< 7 j x (a;)) = 
rlenti(2). Combining this with Lemmas 40 and 41 establishes rlent] (j) = rlenti(ji — 1) for all 
j € {2, 3,..., d — 1}. □ 

We can now use the similarity between A' 2 and the 2 -curve that forms the beginning of any 
approximating curve A' k , to prove the following: 

Lemma 43. rlent^i = rlent 1,1 for all k > 1. 

Proof. For k = 1 the Lemma is trivial. Now consider the case k > 2. By self-similarity, A' k starts 
with a non-reverse isometry of either A' 2 or A! 2 . In the first case we have rlent= rlenti .1 o |< 7 ^i[|; 
in the second case we have rlent k,i = rlexti^ = rlenti p oijj o\a k \\. In either case, rlentfc j i(l,... ,d) 
is a permutation of rlent 11 (1,..., d), which, by Lemma 42 can have only two values: it is either all 


zeros, or rlenti i(d) = 0 and otherwise it is all ones. In the first case, any permutation is without 
effect so rlent ^4 = rlentyi for any k. In the second case, we must have rlentfc j i(d) = 0 for any k 
because e^o = d, and it follows that rlentfe i i(l,..., d — 1 ) is all ones for any k. 

□ 


Lemma 44. The combination of Lemmas \42\ and \^S\ is equivalent to: 

• isneg(f7^" 1 ( t j)) = 0 for all k and all j; or 

• isneg(cr^" 1 (j)) = 0 if k is even or j = d, and isneg(cr A T^(j)) = 1 if k is odd and j < d. 

Proof. Recall Observation 30 rlentfc^j) = (entfcp^+Ufcqfj]) mod 2, where entfc^j) = isneg(cr^ 1 (j)) 
(mod 2). Therefore, isneg(f7^" 1 (j)) = entfcj(j) = rlentfc^j) + Vk,i [j] (mod 2). Since the entry 
point of Ck, 1 is, by definition, v k +i,i, we obtain isneg(cr A T 4 ( ia (j)) = rlentfc + i,i(j) + v k +i,i[j] = 
rlentfc+i,i(j) + isneg((7^^(j)) (mod 2), and therefore sign(t7^ ltl (j)) = sign(f7“j(j)) if and only if 
rlentfc+ipQ') = 0. 

The lemma now follows by straightforward induction from the base case k = 0 (in which 
case <7~f\ is the identity permutation) and the possible values of rlent^i as given by Lemmas 42 
and EH TT 
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Theorem 45. If f is a self-similar, hyperorthogonal, well-folded, space-filling curve mapping [0,1] 
to [ 0 , l] d , then, modulo reflection, reversal and rotation, the entry /( 0 ) is either ( 0 , 0 ) or 

(4 1 0 ) 

Proof. This is a direct translation of Lemma [44] using Observation [ 8 ] □ 


5.6 Construction of self-similar curves 

We will now show that curves with the entry points that may exist according to Theorem |45| do 
indeed exist for any d > 3: 


Theorem 46. For any d> 3, there is a self-similar, hyperorthogonal, well-folded d-dimensional 
space-filling curve starting at ( 0 ,..., 0 , 0 ) and there is a self-similar, hyperorthogonal, well-folded 
d-dimensional space-filling curve starting at ( 3 ,..., 3 , 0 ). 


Proof. It suffices to show that the construction of Theorem 19 


direction (—(d — 1 )), and signs of af\ corresponding to either ( 0 . 
in a self-similar curve. 

Let x be rlenti i i(l). Applying the translation of Lemma 


with entry direction (d), exit 
.., 0 , 0 ) or (|,..., |, 0 ), results 


44 


and Theorem 


45 


in the other 


re relative entry 


direction, we find that both starting points satisfy Lemma [42] and Lemma [43] so t: 
coordinates of the first child curve Ck,i on any level k are given by rlentfcp(j) = x for 1 < j < d — 1 , 
and rlentfc ; i(d) = 0 . 

By Lemma 21 all child curves of the constructed approximating curves A'^ , A\ ,... are an image 
of an isometry of G'(d). By Lemma 36 the relative coordinates of the entries and exits of the 


one-step inflation of each such child curve are permutations of each other, and by Observation 31 


a child curve’s relative entry coordinates are the previous child curve’s relative exit coordinates. 
Thus, the relative entry and exit coordinates of the one-step inflations of all child curves are 
permutations of (x ,..., x, 0). Because Theorem 19 guarantees the continuity of the approximating 
curves, we have rlentfc^j) = 0 if \ek,i-i\ = j and rlentfc^') = x if |e*,^— 11 7 ^ j; similarly, we have 
rlext fc>i (j) = 0 if |e fc)i | = j and rlext fc)i (j) = x if |e fc)i | 7 ^ j. 

Thus, any extended child curve’s inflation, to any depth of recursion, consists of child curves of 
type 0 and 1 , with the entry and exit points determined by the fact that all relative entry and 
exit coordinates are equal to x, except that we have rlent^ (]e*_x|) = 0 and rlext^ (|e^ |) = 0. Thus, 
the entry point (or, in the case of reversal, the exit point) of the inflation of any vertex tq^ to a 
depth of k levels is completely determined by ayy in the same way in which the entry point of Ak 
is determined by 00,1 (which is the identity permutation). As a result, the inflation of iq^ must be 
a translation of o\,i{Ak) or its reverse; hence the space-filling curve is self-similar. □ 


It turns out that there are actually very few such curves for d = 3 and d = 4: 

Observation 47. If d = 3 or d = 4, Lemma |15| leaves no choice with respect to the last two 
elements, the third-last element, and the first element of the permutations |(Tfcy| in a self-similar 
curve. 


Proof. By Lemma |24[ we may assume that the space-filling curve is approximated by extended 
hyperorthogonal well-folded curves A' 0 , AI X ,... with the entry and exit direction fixed at (d) and 
(—(d — 1)), respectively. 

The last two elements of any permutation \<Jk,i\ must be the two different axes of the edges 
incident on Vk,i- 

If d = 3, the third-last (and first) element must be the only remaining axis. 

If d > 3, the third-last element must be the third axis that is within edge distance 1 from 
Vk,i- For i = 1, this third axis is [efc^l, which must differ from le^ol and |efc ; i| = |efc ; 3 1, otherwise 
(efc, 0 ) efe,i, ek, 2 , e^) would constitute a sequence of four edges with only two different axes, con¬ 
tradicting Definition [ 9 ] By a symmetric argument, for i = 2 d * k , the third axis is \ek,i- 2 \- For 
1 < i < 2 d * k , a third axis must also exist, otherwise the two edges preceding Vk,i and the two edges 
following Vk,i would constitute a sequence of four edges with only two axes. If d = 4, with the 
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Figure 5: The three-dimensional, self-similar, hyperorthogonal, well-folded space-filling curve with 
starting points (0,0,0) (a, left) and (|, |,0) (/3, centre), and the three-dimensional curve by Butz 
and Moore (right). The bold grey curve shows A±. The solid black curves depict the child curves 
of Ai, the dashed lines between them indicate how they are connected. The symbols next to 
the child curves indicate whether they are reversed, with arrow, or not, without arrow. For the 
Butz-Moore curve, no such indications are given, because the curve is symmetric and there is 
no need to distinguish between reflections and reversals. The white and black dots indicate the 
location of the entry /(0) and the exit /(1). 

last two elements and the third-last element fixed, the first element must be the only remaining 
axis. □ 


Corollary 48. If d = 3 or d = 4, there are exactly two self-similar, hyperorthogonal, well-folded 
d-dimensional space-filling curves. 


Proof. For self-similar curves, by Lemma |24[ we may assume the entry and exit direction to be 
fixed at (d) and (—(d — 1)), respectively. For the starting point, that is, the signs of v kl (j) for all 
k and j, only two combinations are possible (Theorem |45[). Theorem 19 states that this leads to 
two unique hyperorthogonal, well-folded space-filling curves in which the elements of each |<Xfcy| 
are sorted by order of decreasing local edge distance to Vk,i in A' k . By Observation 47 for d = 3 
and d = 4, there is no other way to order the elements of each |cr fc i |. □ 


The two three-dimensional self-similar, hyperorthogonal, well-folded space-filling curves are 
illustrated in Figure [5j left (a), and centre (/3). 


6 Implementation in software 

6.1 Typical operations 

In order to apply hyperorthogonal well-folded space-filling curves in practical applications, one 
needs to implement one or more operators based on these curves. Recall that the space-filling 
curves under consideration in this paper are functions / : [0,1] —> [0, l] d , with approximating curves 
Aq, Ai,.... Common operators for such curves include: 

• discrete index-to-point conversion: given a resolution parameter k and an index i £ 
{1,..., 2 d * k }, compute the coordinates of vertex Vi of A k ; 

• continuous index-to-point conversion: given a number x £ [0,1], calculate f(x); 

• discrete point-to-index conversion : given a resolution parameter k and the coordinates of a 
vertex v of A k , compute the index i such that Vi = w, 

• continuous point-to-index conversion: given a point p £ [0,1] d , calculate / _1 (p); 
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• discrete comparison: given a resolution parameter k and the coordinates of two vertices u, v 
of A*,, compute which of the two appears earlier along A&; 

• continuous comparison: given the coordinates of two points p, q £ [0, l] d , compute which 
of the two appears earlier along the curve /, that is, decide whether / _1 (p) < f~ 1 (q) or 

f~ l ip) > / _1 («)- 

There is a catch here: the inverse f^ 1 of a space-filling curve / is not immediately well-defined. 
For a given point p, there may be an approximating curve A such that two (or more) hypercubes 
Hi and Hj 1 corresponding to vertices Vi and Vj on A*,, each have p on their boundary, where 
j — i > 1. That implies that there will be a value x £ [(i — l)/2 d * k , i/2 d * k ] and a different value 
y £ [( j — l)/2 d * k , j/2 d * k ] such that f(x ) = f{y) = p. A common solution to obtain a unique value 
for / _1 (p) is to “err on the far side”: for any level k, assign each point p to the vertex Vi of A*, 
whose corresponding hypercube Hi contains the immediate vicinity of p in the direction away from 
the origin. In other words, we define / -1 (p) as the limit of the elements of {x £ [0,1] | f(x) = p'} 
as p' approaches p in a straight line directed towards the origin. A drawback of this solution is 
that / _1 (p) is undefined when one or more of the coordinates of p are equal to 1. An alternative 
solution could be to define / _1 (p) as the smallest value x such that f(x) = p. 

In the context of this publication, it would go too far to go into the details of the optimal 
implementation of each of the operators mentioned above, with various definitions of /“. For¬ 
tunately, the implementations of these operators share the same global structure: starting from 
the unit hypercube, the operator zooms in onto successively smaller hypercubes until the required 
output can be delivered. In Section |6.2| we sketch briefly how to implement the continuous com¬ 
parison operator with the err-on-the-far-side definition of / _1 for the d-dimensional self-similar 
hyperorthogonal well-folded space-filling curves that underlie Theorem [46] Further details are 
provided in Appendix [A] A good understanding of our implementation should enable the reader to 
implement any of the other operators. 

6.2 Implementation of a comparison operator 

It is relatively easy to implement an efficient comparison operator that decides which of any two 
given points comes first along a d-dimensional, self-similar, hyperorthogonal, well-folded space-filling 
curve. For a fixed choice of space-filling curve /, a recursive implementation would take as input 
two points p,q £ [0, l) d that need to be compared, along with a signed permutation a that specifies 
how the given curve is placed in the unit cube, and the direction of the curve (forward or reversed). 
Let S(p) and S(q) be the subcubes of width 1/2 that contain p and q , respectively. 

If p = q, one point does not precede the other. Otherwise, if S(p) ^ S(q ), one can decide 
immediately which point comes first, based on the relative order of the vertices that represent S(p) 
and S(q) along the approximating 1-curve u(G(d)). Finally, if S(p) = S(q), that is, p and q lie in 
the same subcube of width 1/2, then their relative order can be decided by a recursive call with: 

• the points p and g, scaled and translated according to the transformation that maps S(p) to 
the unit cube; 

• the signed permutation and direction that specifies how the space-filling curve traverses S(p). 

In fact, thanks to the structure of the approximating curve a(G(d)), one can examine the coordinates 
of p and q one by one, from the coordinate in dimension |<r(d)| down to the coordinate in dimension 
|cr(l)|: as soon as a coordinate is found in which the binary representations of the fractional parts 
of p and q differ in the first bit, one can decide which of the two points precedes the other. Only if 
p and q are equal in the first bits of all coordinates, the algorithm needs to go in recursion. 

To be able to make the recursive call, the algorithm needs to determine the permutation to 
use in recursion, that is, the transformation that maps the complete space-filling curve / to the 
section within S(p), modulo scaling and translation. For the curves described by the constructions 
of Lemma [18] and Theorem [46] this is relatively straightforward. To determine the unsigned 
permutation to be used in recursion, we sort the d coordinate axes by decreasing local edge distance 
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S(p). This sorted list of axes can be constructed on the fly in 0(d) time while examining the 
d coordinates of p and q to decide in which subcube they lie. By Lemma [18J the sorted list of 
axes gives us the (unsigned) permutation to use in recursion. The signs of the permutation to use 
in recursion now follow from applying the observations on relative entry points and permutation 
signs calculated in the previous section. For further details and pseudocode of a (non-recursive) 
implementation, see the appendix. 

If the binary representations of the coordinates of p and q consist of k bits per coordinate, and 
we can extract these bits in order of decreasing significance in constant time per bit, then the 
complete comparison operator runs in 0(d * k ) time. 


7 Evaluation 

7.1 Comparing to the Butz-Moore curves 

The generalization of Hilbert’s curve to d dimensions by Butz [4], as implemented by Moore (T5j . 
is a self-similar well-folded curve with starting point in the origin, in which the orientations (and 
therefore, the signs of the inverse permutations) of the child curves of A\ are the same as in 
our hyperorthogonal well-folded curves. Concretely, \ci{d)\ = 1 for i G {l,2 d }, and |cr 2 (d)| = 
max(|ei_i|, | |) for 1 <i < 2 d . However, otherwise the permutations are different: all permutations 

in the Butz-Moore curves are rotations (in the permutation sense of the word), so \ui{j)\ = \<Ji(d)\+j 
(mod d). For a graphical description of the 3-dimensional curve, see Figure [5] (right). 

Theorem 49. The d-dimensional Butz-Moore curve contains subcurves with box-to-curve ratio 
H(2 d / 2 ). 

Proof. Assume d > 3. Then G(d) contains a sequence (1,2,—1,(2 + |_c?/2j), 1) or a sequence 
(1, —2, —1, (2 + [d/ 2J), 1). Hence, for the child curves of Ai, there is an i such that |eq(d)| = 2, 
|ej| = 1, and |cri + i(d)| =2+ [d/2\. Now consider the last 2L d / 2 Jedges of Ci and the first 2L d / 2 J- 1 
edges of 6{ + i. By Lemma[4j each of these two sets of edges has \d/2\ different axes. As a result of 
the rotations |<Ti| and |cr i+ i|, these sets of axes include {3,..., 2 + [cZ/2j} and {3 + \ d/2 \,..., d, 1}, 
respectively, where the latter set reduces to {1} if d < 5. Together these sets constitute at least 
the set {1,..., d} \ {2}. Thus the curve through the last 2L d / 2 J — 1 + 1 vertices of Ci and the first 
2 [_d/ 2 \~i _|_ i vertices of C)+i has bounding box volume at least 2 d_1 , and hence the worst-case 
box-to-curve ratio is at least 2 d_1 /(2 d / 2 + 2) = H(2 d / 2 ). □ 

The worst-case box-to-curve ratio of the Butz-Moore curves is thus in sharp contrast with the 
worst-case box-to-curve ratio of our hyperorthogonal, well-folded curves, which have BCR at most 4 
for any d. For verification we also calculated the actual worst-case BCR values for d £ {2,3,4,5, 6} 
with the software from Sasburg m (Table 0. Further investigations may be done into average 
BCR values over curve sections of a given size, both for the hyperorthogonal and the Butz curves. 

It should be noted, however, that BCR may not be the only relevant measure of bounding-box 
quality. Haverkort and Van Walderveen [7] argued that, at least for d = 2, the size of the boundary 
of a bounding box may be as important as its volume—although volume and boundary size are 
usually correlated. Using Sasburg’s software with a generalization of the worst-case bounding box 
perimeter ratio from Haverkort and Van Walderveen to higher dimensions, we found that by this 
measure, already for d = 3, the self-similar hyperorthogonal well-folded curve with starting point 
(|, |, 0) is better than the Butz curve. 

7.2 Lower bounds 

In this work we study space-filling curves that can be described by a series of approximating 
curves Aq, Ai, ..., A n , where A & is a curve on the fc-cube. Within this context, we restricted our 
search for curves with good worst-case BCR first to face-continuous curves; then, more specifically, 
to well-folded curves; then to hyperorthogonal well-folded curves; and finally to self-similar, 
hyperorthogonal, well-folded curves. We found that if d = 3 or d = 4, there are only two self-similar 
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Table 1: Worst-case box-to-curve ratios for various curves in up to 6 dimensions. 


curve 

d = 2 =3 

= 4 

= 5 

= 6 

> 7 

lower bound face-continuous 


2.00 2.54 

3.15 

3.54 

3.76 

4~16/(2 d +3) 

best claimed non-self-sim. 


2.22 a 2.89 b 





self-sim. hyp. well-fid. /(0) = (0,. 

.,0,0) 

2.40° 3.11 

3.53 

3.76 

3.88 

< 4 

self-sim. hyp. well-fid. /(0) = (|,. 


3.14 

3.67 

3.83 

3.92 

< 4 

lower bound non-face-continuous 


3.00 3.50 

3.75 

3.87 

3.93 

4 - 4 /& 

Butz-Moore 


2.40 c 3.11 

4.74 

7.08 

10.65 

n(2 d / 2 ) 


a dQ-cui've ITT] analysed by H&vW [7|; b Iupiter [5]; c Hilbert’s curve sjj 


hyperorthogonal well-folded space-filling curves. For d = 5 and up, there are many more, as 
Lemma 15 then starts to leave room for swaps among the first elements of the permutations <Jk,i- 
We will now address the question of how much room for further improvement there is within these 
restrictions or if some of these restrictions are dropped. 

For d = 2, Haverkort and Van Walderveen [7j report that the BCR of any section of the well- 
folded, non-self-similar /3f2-curve m is 2.22 in the worst case, and for d = 3, Haverkort jjjjj claims 
a fairly complicated, non-self-similar, face-continuous curve with a worst-case BCR of 2.89. These 
two constructions, which do not easily generalize to higher dimensions, constitute improvements of 
less than 10% with respect to the self-similar hyperorthogonal well-folded curves. 

For larger values of d , no face-continuous curve can be much better than any hyperorthogonal 
well-folded curve, since the first is subject to a lower bound that quickly approaches the upper 
bound of the latter as d grows. The proof is based on the fact that any such curve must contain a 
sequence of at most 2 d ~ 2 + 1 edges that have all axes {1,..., d}. 


Lemma 50. Let X be a k-curve constructed by inflation of a single vertex, and let S be a subcurve 
of X. If vol(S')/vol(V) > 2 d ~ 1 /(2 d — 1), then the bounding box of S is the bounding box of X. 


Proof. The proof goes by induction on increasing values of k. 

For k = 0, we have vol(V) = 1 and S only satisfies vol(S')/vol(A) > 2 d ~ 1 /(2 d — 1) if S = X, 
in which case the bounding box of S is indeed the bounding box of X. 

Now suppose the lemma holds for (k — l)-curves, and consider a series of curves Xq, X-\,, Xk 
where A'o is a single vertex and each curve Xi (i > 0) is constructed by inflating A,;_i. Let 
Vi,... ,vd be the vertices of X\, and let Ri be the (k — l)-curve within Xk that results from 
inflating Uj. 

Let S' be a subcurve of Xk with vol(S)/vol(Xfc) > 2 d ~ 1 /(2 d — 1). Since 2 d ~ 1 /(2 d — 1) > 1/2, 
the curve S consists of, at least, a subcurve Y of a curve R y , a subcurve Z of a curve R z , and the 
complete curves Ri for y < i < z, where z ~ y = 2 d ~ 1 . 

We define vol(R) = note that vol(i?i) = vol(R), regardless of i. We have vol(Y) + 

vol (Z) = vol(S) — vol(-Ri) > 2 d ~ 1 /{2 d — 1) * 2 d * vol (R) — (z — y — 1) * vol (R) = 

(2 d * 2 d ~ l /(2 d — 1) — (2 d_1 — 1)) * vol(i?) = (l + 2 d ~ 1 /(2 d — 1)) * vol(i?). Hence, since vol (Z) < 
vol (R), we have vol(T) > 2 d ~ 1 /(2 d — 1) * vol(i?) and thus, vol(Y)/vo^-R^) > 2 d ~ 1 /(2 d — 1). By a 
symmetric argument, vol {Z)/ vol(JL) > 2 d ~ 1 /(2 d — 1). 

Therefore the bounding box of S must contain at least the complete bounding boxes of the 
(.k — l)-curves Ri for y < i < z. Since z — y = 2 d ~ 1 , the vertices v y ,..., v z cannot all lie within a 
(d — l)-dimensional 1-cube, so their bounding box must be the full unit cube, and the bounding 
box of R y ,..., R z must be the full bounding box of X. O 


Theorem 51. If f is a spac e-filling curve approximated by a series of curves A$,... ,Ak within 
the framework of Section 1.3 then f has a section with BCR at least 4 — 16/(2 d + 3). 


Proof. Consider the approximating curve A\ with vertices V\,... ,vd and edges e\,... ,eo- 1 - Let 
2 be the smallest z such that Ui=i \ e i\ = {!,-■-, d}, and let y be the largest y < z such that 
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U i-y \e%\ = {1,... ,d}. By our choice of 2 , we have \et\ 7 ^ \e z \ for all i < z, and by our choice of 
y , we have || 7 ^ \e y \ for all y < i < z. Hence, all vertices v, t for y < i < z must have the same 
coordinates with respect to dimensions \e y \ and \e z \, and therefore lie within a (d — 2 )-dimensional 
hypercube of volume 2 d_2 , so z — y < 2 d ~ 2 . Note that the bounding box of v y ,, v z+ \ has volume 
2 d . 

For a given k, let Ri be the (fc — l)-curve within Ak that results from inflating uiy. Let S be the 
subcurve of Ak that starts with the last \2 d ( k ~ 1 ' ) *2 d ~ 1 / (2 d — 1 )] vertices of R y and ends with the first 
p 2 rf ( fc —!) * 2 d ~ 1 J(2 d — 1)] vertices of R z+1 . We have vol(S) < 2 d ( k ~ 1 '> * (2 d ~ 2 + 2*2 d ~ 1 /(2 d ~ 1))+ 2 
(the +2 results from rounding up). By Lemma 50, the bounding box of S is the bounding box 
of the curves R y ,..., R z +\, which has volume 2 “* 2 d ^ k ~ 1 \ Hence, the box-to-curve ratio of the 
section of / corresponding to S is at least 2 d /(2 d ~ 2 + 2 d / ( 2 d — 1) +2 1-d ( fc-1 )). The limit for k —> 00 
is 4 — 16/ [2 d + 3). □ 

For the specific case of d = 2, Haverkort and Van Walderveen |7j prove a stronger lower bound 
of 2 . 

Now suppose we drop the restriction to face-continuous curves. More precisely, suppose we have 
a space-filling curve approximated by a sequence of curves on the grid Ho, A±, ..., where we allow 
our curves on the grid to have diagonal edges, that is, we allow any edge (v, w) such that w 7 ^ v 
and \w[j] — v[j]| < 1 for all j £ {1,..., d}. In that case, the lower bound becomes even worse: 

Theorem 52. If there is a k and i such that Vk,i and Vk,i +1 differ in at least two coordinates (in 
other words: if there is a diagonal edge), then f has a section with BCR at least 4 — 4/2 d . 

Proof. Consider the m-curves A' and Y that replace and i’k,i+i in Hfc +m . Let S be the 
subcurve of X with volume [vol(AT) * 2 d ~ 1 /(2 d — 1)], ending at the exit point of A, and let T 
be the subcurve of Y with volume [vol(T) * 2 d ~ 1 /(2 d — 1 )], starting at the entry point of Y. By 
Lemma 50 the concatenation of S, (), and T now has bounding box volume at least 4 * 2 d * m , 
while vol(5) + vol(T) < 2 d * m * 2 d /(2 d — 1) + 2. Hence, the box-to-curve ratio of the corresponding 
section of / is at least 4/(2 d /(2 d - 1) + 2 1 ~ d * m ). The limit for m -)• 00 is 4 - 4/2 d . □ 


7.3 Questions for further research 

Note that, as Table [l] shows, at least for d up to 6 the lower bound of Theorem [52] for curves with 
“diagonal edges” is greater than the worst-case BCR of the best hyperorthogonal, well-folded curves, 
and for higher dimensions the difference between the lower bound and the upper bound is less than 
1%. Therefore, in terms of worst-case BCR, little is to be expected from non-face-continuous curves 
based on inflation of /c-cubes for increasing k. 

The question remains whether there are hyperorthogonal curves that are not well-folded, and if 
so, whether such curves would also have good bounds on the box-to-curve ratio. In other words: is 
well-foldedness really required in Theorem |13|? Regardless, Theorem |51| shows that in any case, 
there is not much room for finding curves with a better worst-case BCR within the framework of 
Section 11.51 

Can we find space-filling curves with a better worst-case BCR outside this framework? Peano’s 
space-filling curve and its obvious generalization to higher dimensions are based on approximating 
curves Ai on grids of 3 d * 1 vertices. For these curves in 2, 3, 4, 5, and 6 dimensions, Sasburg’s 
software m reports a worst-case BCR of 2.00, 3.06, 3.64, 3.87, and 3.96 respectively. This may 
serve as evidence that, also for these curves, four is an asymptotic upper bound on the worst-case 
BCR, regardless of d. Note, however, that in higher dimensions, the BCR of these curves seems to 
be slightly worse than the BCR of our hyperorthogonal well-folded curves. 

Departing from the framework of Section [TT3| even further: would it be possible to find space¬ 
filling curves with a better worst-case BCR that cannot be approximated by Hamiltonian paths on 
hypercubic grids? Or are such curves also subject to an asymptotic lower bound of 4? 

One may also ask what lower bounds could be proven in more restricted settings than that 
of Section 1.3. For example, Alber and Niedermeier |T] provide a framework for the description 
of generalizations of Hilbert curves that are self-similar and, in the terminology of Haverkort [5], 
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order-preserving : Ak+i is the concatenation of 2 d scaled, translated, rotated and/or reflected 
but not reversed copies of A^. From Lemma |34| we know that any approximating curve A 2 of a 
self-similar hyperorthogonal well-folded space-filling curves contains child curves of both types 


(zero and one), that is, it contains both non-reverse and reverse isometries of G'(d). By Lemma 25 


these are really different: no non-reverse isometry of an inflation of G'(d) can visit its vertices in 
the same order as G'{d). So no self-similar hyperorthogonal well-folded space-filling curves exist 
without reversal, and thus we get: 


Corollary 53. No d-dimensional self-similar hyperorthogonal well-folded curve for d > 2 can be 
described within the framework of Alber and Niedermeier w- 

Are the curves that can be described within the framework of Alber and Niedermeier subject 
to an exponential lower bound on the worst-case BCR? 
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A Implementation of a comparison operator 

In this appendix we explain how to implement an efficient comparison operator that decides which 
of any two given points comes first along a d-dimensional self-similar hyperorthogonal well-folded 
space-filling curve. Algorithm [l] gives an implementation for a curve with entry point (0,... ,0), 
assuming d > 3. (For d = 2, one could use any implementation of Hilbert’s curve.) We will briefly 
explain how the algorithms works below. We have also tested the algorithm and verified that it 
correctly orders all grid points along hyperorthogonal, well-folded curves, for all grids of 2 d * k points 
with 3 < d < 6 and 2 < k < 12/d. A truly efficient implementation may call for the use of various 
bit tricks (for example, an array whose elements are 1 and — 1 could be encoded as a single binary 
number); however, in the interest of readability, with our implementation we strive to stay closer to 
the theory of Section [5] and avoid tricks that would hide too much of what is going on conceptually. 


A.l Input and output of the repeat and for loops 

We will first describe the input and output of the repeat and for loops. After that we will explain 
how this functionality is implemented. 

The code takes two points p,q £ [0, l) d that need to be compared. In the for loop (Lines 
[5] to 17) the algorithm tries to decide which of the two points comes first along the curve, 
assuming that the curve is reversed as specified by direction (1 means: forward, not reversed; —1 
means: reversed), and rotated and reflected according to the signed permutation cr specified by 
unsgnedPrm and sgnsInvPrm. Here unsgnedPrm[ 1,..., d\ gives the absolute values of er(l),..., a(d) 
and sgnsInvPrm[l ,..., d] holds the signs of cr~ 1 (l),..., cr _1 (d) (the entries unsgnedPrm[ 0] and 
sgnsInvPrm[0] are sentinels that are used to prevent indexing arrays out of bounds on Lines [6] 
and 17 when i = d). On LineJlJ the direction is initialized to forward and a is initialized to the 
identity permutation. 

If p and q lie in the same subcube H of width 1/2, the for loop ends without returning a result, 
but as a side effect, it will have done the following: 

• p and q are scaled and translated according to the transformation that maps H to the unit 
cube; 

• the signed permutation that specifies how the curve traverses H has been determined and 
stored in unsgnedChldPrm and sgnsInvChldPrm (modulo some small “mistakes”, which will 
be corrected in Lines 18 to 24); 


• the position of H in the order in which the curve traverses the unit cube has been stored in 
sbcubeld (0 for the first subcube; 2 d — 1 for the last subcube). 


The algorithm will then, on Lines 18 to 26 correct the “mistakes” and set up unsgnedPrm, 


sgnsInvPrm and direction for the next iteration, which effectively zooms in on the subcube H that 
contains p and q. If, eventually, p and q cannot be distinguished, the algorithm returns 0. 


A.2 


Deciding in which subcube p and q lie 

1 now describe how the for Iood determines in which 


We will now describe how the for loop determines in which subcube(s) p and q he. For now, the 
reader may ignore the assignments to entrAxs , extAxs , unsgnedChldPrm: these have a role in 
determining the signed permutation that specifies how the curve traverses the common subcube (if 


any) of p and q ; we will get back to that in Section A.3 


Recall that the space-filling curve that fills the unit cube is approximated by a curve <r(G{d)) 
with vertices v\,... ,vd- Each vertex Vj corresponds to a hypercube of width 1/2, and in particular, 
there will be two indices j and k such that Vj and i>k correspond to the hypercubes Hj and Hk 
that contain p and q, respectively. Note that, thanks to our decision to “err on the far side”, for 
any m and for any i £ { 1 ,..., d}, the first bit of the fractional part of coordinate i of any point in 
H m is equal to v m [i]. The main goal of the for loop is to identify whether j < k or j > k, and, if 
j = k, what is their value. 


31 









Algorithm 1: Comparison operator based on the d-dimensional self-similar lryperortlrogonal 
well-folded space-filling curve with entry point (0,..., 0), d > 3. 

Input: Points p = (p[ 1],... ,p[d\) and q = (g[l],..., q[d\) in [0, l) d 

Output: —1, 0, or 1: if 1, p precedes q along the curve; if 0, p = q; if —1, p follows q 

1 direction 4— 1; unsgnedPrm[0 ,..., d] 4— [0,..., d]; sgnsInvPrm[ 0,..., d] 4— [1,..., 1] 

2 repeat 

3 entrAxs 4 - unsgnedPrm{d\; extAxs 4— unsgnedPrm[d — 1] 

4 quartAxs -4- unsgnedPrm[d ]; sbcubeld 4- 0 

i 4- 1 to d do 

axis 4- quart Axs; quart Axs 4- unsgnedPrm[d — i\; sbcubeld 4- 2 • sbcubeld 

// figure out in which half of the cube p and q are: 
p[axis\ 4- 2 ■ p[axis]; plnTheBack 4- [p[cms]J; p[axis] 4- p[axis] mod 1 
q[axis\ 4- 2 • q[axis\; qlnTheBack 4- L<?[aads]J; q[axis\ 4- q[axis] mod 1 
if plnTheBack ^ qlnTheBack then 

// on different sides: return 1 if p comes first; — 1 if q comes first 
return direction • sgnsInvPrm[axis] ■ sign(qInTheBack — plnTheBack) 

// determine sign such that entry point lies on outside: 
sgnsInvChldPrm[axis] 4— 1 — 2 ■ plnTheBack 

if plnTheBack = isneg(sgnsInvPrm[axis]) then 

| unsgnedChldPrm\i — 2] 4— extAxs; extAxs 4- axis // p and q in 1st half 

else 

unsgnedChldPrm[i — 2] 4— entrAxs; entrAxs 4- axis // p and q in 2nd half 

sbcubeld 4- sbcubeld + 1 

sgnsInvPrm[quartAxs] 4 - sgnsInvPrm[quartAxs] 

// fill in last two elements of unsigned permutation: 
is unsgnedChldPrm[d — 1] unsgnedPrm[ 1] 

19 unsgnedChldPrm[d] 4- entrAxs + extAxs — unsgnedPrm[l] // the other axis 

// in first and last subcube it is the other way around: 

20 if sbcubeld 4 {0, 2 d — 1} then swap unsgnedChldPrm[d — 1], unsgnedChldPrm[d] 

// correct first element of permutation in last quarter: 

21 if sbcubeld > 2 d then unsgnedChldPrm[ 1] unsgnedPrm[d] 

II correct entry point to be on inside w.r.t. unsgnedPrm[ 1]: 

22 sgnsInvChldPrm[unsgnedPrm[l}] 4 - sgnsInvChldPrm[unsgnedPrm[T\] 

II correct entry point to be on inside w.r.t. orientation of subcube: 

23 orientation 4- unsgnedChldPrm[d\; if sbcubeld ^ {0,2 d — 1} then 

24 sgnsInvChldPrm[orientation\ 4 - sgnsInvChldPrm[orientation\ 

25 unsgnedPrm = unsgnedChldPrm; sgnsInvPrm = sgnsInvChldPrm 
//if type 1, reverse direction: 

26 if extAxs = orientation then direction 4 - direction 

27 until p = q 

28 return 0 



// p and q are equal 





To this end the for loop implicitly maintains a lower bound lowbnd and an upperbound uppbnd 
on j and k. In successive iterations, the gap between these bounds is narrowed until we either find 
j ^ k , or lowbnd = j = k = uppbnd. In the last case, sbcubeld eventually holds the value of j — 1 
(this is because this article generally indexes vertices starting from one, but the implementation 
starts from zero). Specifically, the following invariant is valid just before each execution of Line [ 6 ] 
lowbnd = 2 d+1 ~ l * sbcubeld + 1 < j < k < 2 d+l ~ l * (sbcubeld + 1 ) = uppbnd. Just after LinejbJ the 
following holds: (i) lowbnd = 2 d ~ l * sbcubeld + 1 < j < k < 2 d ~ l * (sbcubeld + 2 ) = uppbnd ; (ii) 
axis = | cr(d + 1 — z)|, and (iii) quartAxs = \a(d — z)|. 

Note that due to the properties of G(d), we always have that vi ow b n d > • •., v up pbnd is a translation 
of cr(G(d + 1 — i)) or its reverse, and this curve consists of the concatenation of a(G(d — i)), an 
edge with axis a(d + 1 — *), and a(G(d — i)). In iteration i of the loop, the algorithm decides 
whether j and k lie in the first or in the second half, that is, before or after the edge with axis 
axis = | cr(d + 1 — i)|. Line [7] reads and removes the first bit of the fractional part of p[axis] (that 
is, Vj[axis]), shifting the remaining bits left for the next iteration of the repeat loop. The bit that 
is read is stored in plnTheBack. Similarly, Line [8] reads and removes the first bit Vk[axis] of the 
fractional part of q[axis\. 

If p and q differ in the bits just read, we can now decide, on Line |10[ which of the two comes 
first along the space-filling curve. In the absence of any reflections or reversals, we would return 1 
if p has the smaller coordinate and — 1 if <7 has the smaller coordinate. However, if this portion of 
the curve is reflected in this coordinate, or if it is reversed, the return value is modified accordingly 
by multiplying with sgnsInvPrm[axis\ and direction. 

If, on the other hand, p and q have the same initial bit in dimension axis , their shared 
coordinate is effectively stored in sgnsInvChldPrm : if p and q lie in the back (the bits read were 
ones), we store —1; if p and q lie in the front (the bits read were zeros), we store 1. Thus, 
isneg (sgnsInvChldPrm[axis]) = plnTheBack = qlnTheBack. The next iteration of the for loop 
must now zoom in onto the isometric copy of a(G(d — i)) that contains Vj and Vk out of the two 
copies that appear before and after the edge with axis axis = \a(d + 1 — i)| in cr(G(d + 1 — i)). 
There are two possibilities: 


Vj and Vk lie in the first part (traversed according to a(G(d — *))): in the absence of reflections, 
this is the case if Vk[axis\ = 0, but if a encodes a reflection in dimension axis, then Vj and Vk 
lie in the first part if Vk[axis] = 1. Line 12 checks for this: the condition evaluates to true 
if and only if Vj and Vk lie in the first part. If so, the upper bound on j and k needs to be 
lowered. This is realized by incrementing the loop counter i and doubling sbcubeld on Line [ 6 ] 
of the next iteration. 


Vj and Vk lie in the second part, traversed according to a(G(d — i)), or equivalently, <r(G(d—i)) 


reflected in coordinate a(d — i). Line 17 implements that reflection, or rather, it “falsifies” 
sgnsInvPrm[a(d—i)\ such that it tricks the next iteration of the for loop into acting according 
to a reflection in coordinate a(d — i). Additionally, the lower bound on j and k needs to be 
raised. This is realized by incrementing sbcubeld on Line |16[ followed by increasing the loop 
counter i and doubling sbcubeld on Line [ 6 ] of the next iteration. 


Note that, when i becomes d+ 1 so that the for loop terminates, as a result of the loop invariant, 
if p and q have the same initial bit in each dimension, sbcubeld will hold the correct value of j — 1 
and fc — 1 once the for loop terminates. 


A.3 Sorting axes by local edge distance 

If the for loop completes, that is, j = k , we need to set up the permutation a to use in the next 
iteration of the while loop. In what follows, we will continue to use cr for the permutation used in 
the current iteration, and we will use o' to denote the permutation to use in the next iteration. 
We will use the construction of Section 14.21 the absolute values of the elements of a' are sorted in 
order of decreasing local edge distance to Vj. To realize this, the algorithm exploits the following 
property of gray codes: 
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Observation 54. Let C be the concatenation of (a), G(m) and (z), where |a| > to, \z\ > m and 
\a\ 7 ^ \z\. If v is a vertex in the first half of G(m), then ed(C,v,z) > ed(C,v,e) for any edge (e) 
in the concatenation of (a) and G(m). Symmetrically, if v is a vertex in the second half of G(m), 
then ed (C,v,a) > ed(C, v,e) for any edge (e) in the concatenation of G{m) and (z). 


Recall from Section |5.2| that a self-similar hyperorthogonal well-folded curve must have an 
extended approximating curve A' = a(A , 1 ) with entry edge (cr(d)} and exit edge (a(—(d— 1)). As 
explained above, in iteration i of the for loop, we are trying to locate Vj and Vk in a section of 
A! that is isometric to cr(G(d + 1 — *)), and we decide whether Vj and Vk appear in the first or 
in the second half of that curve. As a result, in each iteration of the for loop we may be able to 
apply Observation |54| to determine one more axis in the sequence of axes sorted by decreasing 
local edge distance to v :l in A'. Considering this idea more carefully, we see that in the first two 
iterations Observation 54 cannot be applied since some or all of the preconditions |a| > to, \z\ > m, 
and |a| 7 ^ \z\ are violated; this is consistent with the fact that after the last iteration, two axes 
must remain that both have local edge distance zero and cannot be sorted. In the third iteration, 
Observation 54 can be applied if Vj is in the first three quarters of A' , but things go wrong if 


Vj is in the last quarter of A! , which is a reflection of a{G(d — 2)) preceded and followed by an 
edge (a(—(d— 1 ))) (however, in that case, it is clear that |cr(d)| is the axis with the largest edge 
distance). In each of the iterations after the third we can always determine one more axis in the 
sorted sequence. 

In our implementation, the sorting is implemented by assignments to unsgnedChldPrm , sup¬ 
ported by assignments to entrAxs and extAxs. For ease of implementation, an axis is assigned to 
unsgnedChldPrm[i — 2] in each iteration i of the for loop (on Line 13 or 15), but the assignments 


in the first two iterations (to unsgnedChldPrm[— 1] and unsgnedChldPrm[0]) are meaningless and 
without consequence. Throughout the iterations of the for loop, the algorithm keeps track of the 
axes entrAxs and extAxs of the edges that precede and follow the curve cr(G(d + 1 — i)) currently 
under consideration, by the assignments on Lines nnn and |15| Thus, when the for loop ends, 
entrAxs and extAxs store the two axes at edge distance zero to Vj. As noted above, if Vj is in 
the last quarter a wrong assignment to unsgnedChldPrm[ 1] is made in the third iteration; this is 
corrected on Linel 2 Tl 

To complete the permutation a' to use in recursion (modulo the signs), we need to fill in 

unsgnedChldPrm[d — 1] and unsgnedChldPrm[d\ with the axes |cr(e*_ 1 ) | and |cr(ej)|, one of which 

equals cr(l) (by Lemma [ 2 ]). Since the space-filling curve has entry point (0, ...,0), we have 


rlenti(i) = 0 for all i £ {1,..., d}. Hence, by Lemmas 33 to 35. \o'(d — 1)| = er(l) and |er'(d)| is 
the other axis out of |cr(e^—i)| and \a(d)\, unless j £ { 172 “} (and thus, sbcubeld £ {0, 2 d — 1}), in 
which case |cr'(cZ)| = <r(l) and \a'(d— 1)| is the other axis. Corresponding assignments are made in 
Lines [18] to [20[ 


A.4 Computing the signs of the permutation a' of the common subcube 

of p and q 

The correctness of the assignments to sgnsInvChldPrm can be verified by developing the calculations 


of Lemmas 33 to 35 further. Some waypoints for these calculations are the following. Let i be the 
rank of the subcube that contains p and q along the curve, that is, i = sbcubeld + 1 . As before, we 
will continue to use a for the permutation used in the current iteration, and we will use a' = a o cq 
to denote the permutation to use in the next iteration. The assignments on Line m set the signs of 
a 1 such that the relative coordinates of the entry point of C,;, as defined in Section |5~3[ are all zero 
(recall the re lati on between subcube coordinates, relative entry coordinates and signs as given by 


Observation 
and 
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flip t 
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which implies rlentj(j) = 0 if and only if isneg (o i (j)) = Vi[j] mod 2). Lines 


22 


re relative entry coordinates in dimension <r(l) and (if 1 < i < 2 d ) dimension |er'(d]"[7 


Finally, Line [26] reverts the direction in all subcubes of type 1, which means that for subcubes of 
type 1 , the aforementioned settings of the relative coordinates are actually for the exit point, not 
the entry point. We will now explain how the correctness of the settings for subcubes of type 0 can 


be derived from Lemmas 33 to p5] the correctness of the settings of the relative exit coordinates for 
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subcubes of type 1 can be derived similarly, using the relations between entry and exit coordinates 
from Observation [3l1 

Note that for the first subcube, all relative entry coordinates are zero. From Lemmas |33| to [35 
we get that Ci has type 0 if and only if i = 2 d , or if i is odd and i / 1. The correct settings of the 
relative entry point coordinates for the case i = 2 d can now be verified directly using Lemma 35 


For the case of odd i ^ 1, the relative entry coordinates are the relative exit coordinates for Cj_ 1 , 
which differ from the relative entry coordinates for Cj_i only in the orientation |er(cq_i(<i))|, which, 
by Lemma 34 equals |<r(cq(d))| = |<r'(<i)| and differs from er(l). It follows from Lemma 34 that the 
relative entry coordinates for Ci are as follows: rlenti(cr(l)) = rlent 2 _i(cr(l)) = 1; for 1 < j < d 
and j ^ \a '(d)| we have rlenti(cr(j)) = rlenti_i(er(j)) = rlentficfij)) = 0; for j = d and j ^ |er'(d)| 
we have rlentj(cr(d)) = rlenti_i(cr(d)) = rlenti(cx(l)) = 0; and for j = |cr , (d)| (and hence, j ^ 1), 
we have rlenti(a(j)) = 1 — rlenti_i(cr(j)) = 1. This is exactly what the algorithm establishes. 


A.5 Running time 

If the binary representations of the coordinates of p and q consist of k bits per coordinate, then 
the algorithm runs in 0(d ■ k) time (that is, linear in the input size), provided that the necessary 
operations to extract a single bit from a coordinate (Lines [T] and [8]) run in constant time. 


A.6 Variant for a space-filling curve with entry point on face 



rlenti(d) = 0. 
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